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1.0 


Atmospheric  Density  Models 


Knowledge  of  neutral  upper  atmospheric  densities  is  an  important  requirement 
for  the  understanding  and  analysis  of  many  phenomena  under  study  at  AF6L. 
Researchers  in  many  areas  desire  best  estimates  of  composition  and  temp¬ 
eratures  as  inputs  to  models  and  data  analysis  programs.  Prominent  examples 
are  analyses  of  auroral,  airglow  and  ionospheric  measurements. 

The  continuing  need  for  accuracy  in  satellite  tracking  and  ephemerides  predic¬ 
tion  also  results  in  the  need  for  improved  modeling  of  thermospheric  density 
variations.  A  case  in  point  concerns  the  SKYLAB,  whose  orbit  decayed  at  a 
faster  rate  than  predicted,  due  to  the  substantial  increase  in  solar  activity 
in  the  current  cycle.  Recently  the  solar  flux  reached  such  high  levels  as  to 
cause  temporary  loss  of  tracking  of  an  800km  satellite.  Rapid  density  varia¬ 
tions  during  magnetic  storms  continue  to  limit  satellite  tracking  and  predic¬ 
tion  accuracy,  particularly  at  lower  altitudes. 

Density  modeling  research  has  thus  continued  at  a  high  level.  During  the 
period  of  this  contract,  updates  were  made  to  the  Jacchia  19771  and  MS  I 
models,  based  on  continued  analysis  of  mass  spectrometer  data.  These  will 
undoubtedly  be  further  updated  as  new  data  becomes  available.  In  particular 
discrepancies  have  been  found  between  MSIS  temperatures  and  recent  measure¬ 
ments  3  which  may  be  reduced  by  the  use  of  revised  and  new  incoherent  scatter 
data  from  the  Millstone  Hill  Radio  Observatory. 

Support  of  AFGL  efforts  has  included  the  following: 

1.  Development  of  a  software  package  of  existing  recent  models  for 
convenient  usage  by  researchers  in  both  density  modeling  and  other 
fields. 

2.  Development  of  analytic  and  condensed  versions  of  the  Jacchia  1977 
density  model . 

3.  Comparative  evaluations  of  satellite  positional  prediction 
accuracies . 

4.  Studies  of  density  variations  by  analysis  of  Doppler  beacon  tracking 
data. 
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1.1 


Density  Model  Package 

Program  DENMOD  provides  a  convenient  means  to  compute  neutral  thermospheric 
densities  and  temperatures  from  recent  models.  Input  may  be  from  cards  or  a 
satellite  ephemeris  file.  Output  consists  of  a  printout  with  options  for 
storage  on  other  devices. 

As  shown  in  Figure  1,  the  package  consists  of  the  main  program  DENMOD,  an 
output  routine  OUTDEN  and  the  density  computation  software  headed  by  sub¬ 
routine  DENS.  The  latter  software  may  by  itself  be  easily  incorporated  into 
other  programs  if  users  desire.  The  appropriate  calling  parameters  are 
defined  in  the  listings.  The  functions  of  the  subroutines  shown  in  Figure  1 
are  described  in  Section  1.1.2. 
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1.1.1  Functional  Description 

Program  DENMOD  allows  one  to  compute,  for  selected  time  and  positions, 
any  or  all  of  the  following  properties 

Mass  Density 

Constituent  densities 

Exospheric  temperature 

Local  temperature 

Solar  and  geophysical  parameters 

The  following  models  are  available: 

1.  U.S.S.R.  -  Cosmos^ 

2.  Jacchia  19771 

3.  MSIS  longitude  averaged? 

4.  MSIS  with  longitude/UT  variations^ 

The  user  provides  time  and  position  input  in  either  of  two  modes. 

1.  Ephemeris  file  (TAPE1 )  in  LOKANGL  format^  (refer  also  to 
section  2.1.5) 

2.  Cards.  In  this  case  a  range  of  heights  may  be  specified 
for  a  single  time  on  one  input  card. 

In  either  mode  the  user  has  the  option  to  directly  input  geophysical 
information  (solar  and  geomagnetic  activity)  or  have  the  program  obtain 
this  from  tables  on  TAPE2. 

In  addition  to  printouts,  the  results  optionally  may  be  saved  on  binary 
(TAPE3)  and/or  card  image  BCD  (TAPE4),  the  latter  being  convenient  for 
transmittal  to  non-CDC  installations. 

For  the  Jacchia  1977  model  a  special  mass  storage  input  file,  JSDM,  is  required, 
which  contains  tables  of  temperatures  and  constituent  densities  an  functions  of 
altitude  and  exospheric  temperature. 

Figure  2  and  Table  1  define  the  inputs  required  and  the  output  files  produced. 

A  sample  printout  is  shown  in  Figure  3. 
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PROGRAM  OENMOO( INPUT .OUTPUT , TAPE  1 , TAPE2 . TAPE3 , TAPE4 , JSDM) 

C 

C  LOG l CON, INC.  AUGUST, 1979 

C  DENSITY  MODEL  PACKAGE  FOR  FOLLOWING  MODELS 

C  USSR-COSMOS  (ElYASBERG.  ET.  AL. ,  SPACE  RES.  XII,  P.  727,  1972) 

C  JACCHIA  1977  (JACCHIA,  SAO  SPECIAL  REPORT  #375,  1977) 

C  MS  IS  LONGITUDE  AVERAGED  (HEDIN,  ET.  AL.,  J.  GEOPHYS.  RES.  82, 

C  P.  2139,  1977) 

C  MS IS  WITH  LONGITUDE/UT  VARIATIONS  (HEDIN,  ET.  AL. ,  J.  GEOPHYS.  RES. 

C  84 ,  P.  1,  1978) 

C 

C  TAPE1  ■  EPHEMER l S  FILE 

C  TAPE2  ■  AP/F 10.7  FILE 

C  TAPE3  «  BINARY  OUTPUT  FILE 

C  TAPE4  «  BCD  OUTPUT  FILE 

C  JSDM  «  JACCHIA  1977  MASS  STORAGE  FILE 

C 

DIMENSION  aPF10(3) , I ONPAR ( 2 ) 

DIMENSION  I0L(5,2),TLIM(2),A(253),IA(253) 

EQUIVALENCE! A , I A ) 

C0MM0N/0UT0EN/ID(5) ,SEC,P0S{3) , DEN ( 1 5 ,4) , JOUT( 1 1 ) , JM0D(4) , TITLE (8) 
1 .NLINE.UOUTP, JIN 
DIMENSION  L I  ST ( 6 ) 

DATA  LIST/3«-1 ,0,2  —  1/ 

DATA  TLIM/O. , 1 .ElO/ 


INPUT  PARAMETERS 
CARO#  COL 
1  1-80 

2  3 

2  6-16 


2  18-20 

2  28-31 


FORMAT  VARIABLES 
8A1 0  T I T  LE( 1 -8 ) 

PRINTOUT  PAGE  TITLE 
II  JIN  TIME/POSITION  INPUT  MODE 

1 -EPHEMER I S  FILE  (TAPEl),  2-CARDS 
1111  JOUT (1-11) 

OUTPUT  INDICATOR  FOR  ITH  VARIABLE 
1 » YES , 0 -NO 

1.  MASS  0ENSITY(GM/CM3) 

2.  02  DENSl TY ( /CM3 ) 

3.  0  DENSITY  (/CM3) 

4.  N2  DENSl TY(/CM3) 

5.  N  OENSiTY  (/CM3) 

6.  HE  OENSl T Y ( / CM3 ) 

7.  AR  DENSl TY ( /CM3 ) 

8.  H  DENSITY  (/CM3) 

9.  EXOSPHERIC  TEMPERATURE(DEG  K)' 

10.  LOCAL  TEMPERATURE  ( pEG  K) 

11.  GEOPHYSICAL  DATA  AND  SOLAR  EPHEMERIS 
GEOMAGNETIC  INDEX(AP) 

INSTANTANEOUS  SOLAR  FLUX(F10.7) 
SMOOTHED  SOLAR  FLUX(F10.7  BAR) 

SOLAR  RIGHT  ASCENSION( DEG) 

SOLAR  DECLINATION(OEG) 

13  JOUTP 

OTHER  OUTPUT  FILES 

o-none,i=bcd  only, 2«binary  Only , 3-both 
411  JM0D(1-3) 

MODEL  INDICAT0RS(1«YES,0*N0) 

1.  USSR-COSMOS 


Figure  2a.  OENMOD  Input/Output 
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uuuuuuuouuuuuuuuuuuouuuuuuuuuuuuyuuuauuouooouuuuuuuuuuuuo 


2.  JACCHIA  1977 

3.  MS  IS  LONGITUDE  AVERAGED 

4.  MS  IS  WITH  L0NG1 TUDE/UT  VARIATIONS 

(FOR  JIN* 1 ) 


3 

1-3 

13 

IDL{ 1 ,1) 

MONTH 

EARLIEST  TIME 

3 

4-6 

13 

10L( 2 , 1 ) 

DAY 

TO  COMPUTE 

3 

7-11 

15 

I D  L  (  3 , 1  ) 

YEAR 

DENSITIES 

3 

12-14 

13 

I DL ( 4 , 1 ) 

HOUR 

3 

15-17 

13 

IDL(5, 1 ) 

MINUTE 

IF  IDL(I.I)  .LE.  0  COMPUTATION 
STARTS  WITH  1ST  DATA  RECORD 
ON  EPHEMERIS  FILE. 

3 

21-23 

13 

IDL(1.2) 

MONTH 

LATEST  TIME 

3 

24-26 

13 

IDL( 2 , 2 ) 

DAY 

TO  COMPUTE 

3 

27-31 

15 

IDL(3 ,2 ) 

YEAR 

DENSITIES 

3 

32-34 

13 

1  DL(  4 , 2  ) 

HOUR 

3 

35-37 

13 

IDL( 5,2) 

MINUTE 

IF  IOL( 1 ,2)  .LE.  0  COMPUTATION 
ENDS  WITH  LAST  DATA 
RECORD  ON  FILE 


3  41-70  3F10.3  APF10(1-3) 

1.  INSTANTANEOUS  SOLAR  FLUX(F10.7) 

2.  SMOOTHED  SOLAR  FLUX(F10.7  BAR) 

3.  AP 

IF  APF10(1)  .LE.  0.  ALL  3  VALUES  WILL  BE 
OBTAINED  FROM  TABLES  ON  TAPE2. 


(FOR 

3-N 

JIN-2) 

2-3 

12 

IMO(  -t  TERMINATES) 

3-N 

5-6 

12 

IDA 

3-N 

8-11 

14 

IYR 

3-N 

13-M 

12 

IHR 

3-N 

16-17 

12 

IMN 

3-N 

18-23 

F6. 1 

P0S(1)  -  GEOCENTRIC  LAT I TUDE ( DEC  N) 

3-N 

24-29 

F6.1 

P0S(2 )  -WEST  LONGITUDE(DEG) 

3-N 

30-36 

F7.1 

POS ( 3 )  -  LOWEST  HEICHT(KM) 

3-N 

38-42 

F5.3 

OH  -  HEIGHT  INCREMENT(KM) 

3-N 

43-47 

15 

NH  -  NUMBER  OF  HEIGHTS 

3-N 

49-53 

F5. 1 

APF 1 0 ( 1 )  -  INSTANTANEOUS  SOLAR  FLUX(F10.7) 

3-N 

55-59 

F5.1 

APF 1 0 ( 2 )  -  SMOOTHED  SOLAR  FLUX  (F10.7  BAR) 

3-N 

61-65 

F5.1 

APF 1 0 ( 3 )  -  AP 

3-N 

66-70 

15 

IF  APFtO(l)  .LE.  0.,  THEN  ALL  3  VALUES 
WILL  BE  TAKEN  FROM  TAPE2  TABLES. 

IPG  START  OUTPUT  ON  NEW  PAGE  IF  »1 

TAPE  1  INPUT— SEE  AFGL-TR-78-0204  ( LOGICON , INC FINAL  REPORT , 1 978) , P19 


TAPE2  INPUT— 1  RECORD  PER  DAY  STARTING  AT  LEAST  85  DAYS  BEFORE 
EARLIEST  DENSITY  COMPUTATION  TIME 


WORD  NUMBER 

COLUMNS 

FORMAT 

DESCRIPTION 

1 

1-8 

F8,0 

AP  FOR  0-3  HRS  UT 

2 

9-16 

F8.0 

AP  FOR  3-6  HRS 

3 

17-24 

F8.0 

AP  FOR  6-9  HRS 

4 

25-32 

F8.0 

AP  FOR  9-12  HRS 

5 

33-40 

Ps.o 

AP  FOR  12-15  HRS 

Figure  2b.  DENMOD  Input/Output  (Continued) 


r 


C 

6 

41-48 

F8.0 

AP  FOR  15-18  HRS 

c 

7 

49-56 

F8.0 

AP  FOR  18-21  HRS 

C 

8 

57-64 

F8.0 

AP  FOR  21-24  HRS 

c 

9 

65-72 

F8.0 

DAILY  F10.7  AT  17  HRS 

C 

c 

10 

73-77 

15 

MODIFIED  JULIAN  DAY 
( 1/1/73*41683) 

JSDM  INPUT — SEE  DOCUMENTATION 


TAPES  BINARY  OUTPUT 

RECORD  t 

DATA  RECORDS 

WORD  NUMBER 
1 
2 

3 

4 

5 

6 
7 
B 
9 

10 

11  THRU  10+NIMOD 


1 1+NIMOD 

1 i +n IMOD+ i  thru 
1 1+NIMOD+NJOEN 

1 1+NIM00*NJDEN+1 
THRU 

1 1+NIM0D+NJDEN+NJDEN*NIM0D 


HEADER:  (TITLE(I),I«1,8) 


DESCRIPTION 
MQNTH( integer) 

DA  Y ( INTEGER) 

YEAR( INTEGER  -  LAST  2  DIGITS) 

HOUR ( INTEGER ) 

MINUTE! INTEGER) 

SECONDS! REAL ) 

GEOC.  LATITUDE. DEG  N( REAL) 

LONGITUDE, DEG  W( REAL ) 

HEIGHT ,KM( REAL) 

NIMOD  -  NUMBER  OF  MODELS  IN  OUTPUT 

(  INTEGER) 

(IMOD! I).I*1 .NIMOD)  -  CODES  IDENT¬ 
IFYING  THESE  MODELS  (INTEGER)  (SE^ 
BELOW) 

NUDEN  -  NUMBER  OF  VARIABLES  IN  OUTPUT 

(  INTEGER) 

(UDEN! J) ,U*1 ,NJDEN)  -  CODES 
IDENTIFYING  THESE  VARIABLES 
(INTEGER)  (SEE  BELOW) 

( { A ( UD£N ( U ) , IM0D(I ) ) . J>1 , NJDEN ) , 

1.1, NIMOD)  -  RESULTS  FOR  THESE 
VARIABLES  AND  MODELS(REAL) 


TAPE4  BCD  OUTPUT 

RECORD  1  SAME  AS  FOR  BINARY ( 8A10 ) 

DATA  RECORDS 

SAME  AS  FOR  BINARY.  EACH  SPLIT  INTO  3  OR  MORE  LINES 


LINE  NUMBER 
1 

2 

3  THRU 

2+(NIM00»NJ0EN)/7 


WORDS 

1-9 

io  thru  i i+nimod+nuden  : 

1 2+N IMOO+NUDEN 
THRU 

1 1+NlMOO+NjDEN+NjDEN^NlMOD 


format 

12X, 2 ( 12, 1 H/ ) , 12.313, 
F7.1 ,F6.1 ,FB.1 
21  13 

1P7E1  1  .3 


MODEL  codes 

1 .  USSR-COSMOS 

2.  JACCHIA  1977 

3.  MS  I S  LONGIT.  AVERAGED 

4.  MSIS+LONGI TUOE/UT 


VARIABLE  CODES 

1 .  MASS  DENSITY 

2.  02  DENSITY 

3.  0  DENSITY 

4.  N2  DENSITY 

5.  N  DENSITY 

6.  HE  DENSITY 

7.  AR  DENSITY 

8.  H  DENSITY 

9.  EXOSPHERIC  TEMPERATURE 

10.  LOCAL  TEMPERATURE 

11.  GEOMAGNETIC  INDEX 

12.  INSTANTANEOUS  SOLAR  FLUX 

13.  SMOOTHED  SOLAR  FLUX 

14.  SOLAR  RIGHT  ASCENSION 

15.  SOLAR  DECLINATION 


Figure  2c.  DENMOD  Input/Output  (Continued) 
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Record  Number  1 


Word  # 

Symbol 

Description 

1 

LABEL 

BCD  label 

2-9 

TITLE 

80  Character  BCD 

10 

TMINMOD 

Minimum  exospheric 
temp,  in  table 

11 

TMAXMOD 

Maximum  exospheric 
temp. 

12 

NTSTEP 

Number  of  Temperature 
steps 

13 

UTMOO 

=  (TMAXMOD  -TMINMOD) 
NSTEP 

Temperature  step  size 

14 

NR EG  ( <7  ) 

Number  of  height  regions 
(=number  of  subsequent 
records) 

15 

NSUM 

Unused 

16-23 

ALZSTEPi 

Natural  logs  of  boundaries 
(km)  of  height  regions. 

24-30 

NZSTEP-j 

Number  of  height  steps  in 
each  region 

31-37 

ALDZ-j 

Step  size  in  natural  log 
of  height  in  each  region. 

(ALZSTEPm-ALZSTEPi) 


NZSTEP-j 


38 

GNMOD 

Gravitational  acceleration 
at  Earth  surface  (km/sec2) 

39 

RNMOD 

Earth  radius  (km) 

40 

RSTAMOD 

Universal  gas  constant 

41 

AVOGMOD 

Avogadro's  number 

42 

NSPMOD  (=7) 

Number  of  species 

43-49 

AMWMOD-j 

Molecular  masses  of 
species 

50-56 

MODRLEN j 

]= (NSPMOD+1 ) 

Lengths  of  subsequent 
records 

*(NTSTEP+1 ) 
*(NZSTEPj+l )] 

The  species  will  be  given  in  the  order: 
0^,  0,  N2,  N,  He  Ar,  H 


Table  1.  Format  of  Hass  Storage  File  "JSDM" 
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Record  Number  IREG+1  1<IREG<NREG) 


Word  # 


Description 


1 

2  thru 
NSPMOD+1 


Local  temperature  of  height  exp  [ALZSTEP(IREG)] 
and  exospheric  temperature  TM1NM0D 
Logs  base  10  of  number  densities  (m”3)  for  molecular 
species  at  this  height  and  exospheric  temperature. 


NSPMOD+2 

thru 

(HTSTEP+1 )* 
(NSPMOD+1 ) 


Repeat  of  words  1  thru  NSPMOD+1  for  this  height  and 
equally  spaced  exospheric  temperatures  TMINNOD  +  DTMOD 
thru  TMAXMOD 


(NTSTEP+1  )*(NSPM0D+1 )  +  l 
thru 

MODRLEN  (I REG) 


Repeat  of  words  1  thru  (NTSTEP+1 )*(NSPMQD+1 )  for 
remaining  heights 
exp  [ALZSTEP  (IREG) 

+  I  *ALDZ  (IREG)] 

I  =  1  thru  NZSTEP  (IREG) 


Table  1.  Format  of  Mass  Storage  File  "JSDM"  (continued) 
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EXCEDE- II  SPECTRAL 


MOOElS:  2= J ACCH I A  1977  3=MSIS  LONG  AVERAGEO 


MO/DA/YR 

HR 

MN 

SC 

LAT 

W.  LON 

HEIGHT 

MASS  DEN  02  DEN 

0  DEN 

N2  DEN 

AR  DEN 

TLOC 

MODEL 

DEG 

DEG 

KM 

GM/CM3  /CM3 

/CM3 

/CM3 

/CM3 

OG  k 

10/19/79 

5 

47 

0 

G5. 1 

147.5 

107.0 

1 .7916-10  5 . 672E+1 1 

2. 996E+1 1 

2 . 996E+1 2 

2 . 520E  +  1 0 

226 

2 

10/19/79 

5 

47 

0 

65.1 

147.5 

107-0 

HEIGHT  OUT  OF  RANGE 

FOR  MODEL 

3 

10/19/79 

S 

47 

0 

65.1 

147.5 

108.0 

1  .  505E-I0  4 . 589E+1 1 

2. 791 E+1 1 

2 . 523E+1 2 

1 . 999E+ 1 0 

233 

2 

10/19/79 

5 

47 

0 

65.1 

147.5 

108.0 

HEIGHT  OUT  OF  RANGE 

FOR  MODEL 

3 

10/19/79 

5 

47 

0 

65.1 

147.5 

109.0 

1 • 266E-10  3 . 708E+1 1 

2.590E+1 1 

2. 1 28E+12 

1 . 592E+1 0 

242 

2 

10/19/79 

5 

47 

0 

65.1 

147.5 

109.0 

HEIGHT  OUT  OF  RANGE 

FOR  MODE  L 

3 

10/19/79 

5 

47 

0 

65.1 

147.5 

110.0 

1 . 0G8E- 1 0  2 . 995E+ 1 1 

2.395E+1 1 

1 . 799E+1 2 

1  . 273E+1 0 

251 

2 

10/19/79 

5 

47 

0 

65 . 1 

147.5 

110.0 

HEIGHT  OUT  OF  RANC-E 

FOR  MODEL 

3 

10/19/79 

5 

47 

0 

65.1 

147.5 

111.0 

9.039E-11  2.425E+11 

2 . 206E+1 1 

1 . 526E+1 2 

1 . 024E+1 0 

262 

2 

10/19/79 

5 

47 

0 

65.1 

147.5 

111.0 

HEIGHT  OUT  OF  RANGE 

FOR  MODEL 

3 

10/19/79 

5 

47 

0 

65.1 

147.5 

112-0 

7.682E-11  1.972E+11 

2 . 027E+1 1 

1 . 299E+1 2 

8 . 282E+09 

273 

2 

10/19/79 

5 

47 

0 

65.1 

147.5 

112.0 

HEIGHT  OUT  OF  RANGE 

FOR  MODEL 

3 

10/19/79 

5 

47 

0 

65.1 

147.5 

113-0 

6.S61E-11  1.614E+11 

1 .8S8E+1 1 

1.11 OE+1 2 

6.743E+09 

285 

2 

10/19/79 

5 

47 

0 

65.1 

147.5 

113.0 

HEIGHT  OUT  CF  RANGE 

FOR  MODEL 

3 

10/19/79 

5 

47 

0 

65.1 

147.5 

114.0 

5.635E-11  1.333E+11 

1 . 702E+1 1 

9.539E+1 1 

5 . 529E+09 

297 

2 

10/19/79 

5 

47 

0 

65.1 

147.5 

114-0 

HEIGHT  OUT  OF  RANGE 

FOR  MODEL 

3 

10/19/79 

5 

47 

0 

65.1 

147.5 

115.0 

4.868E-1 1  1  -  1 12E  +  1 1 

1 . 559E  +  1  1 

8.238E+1 1 

4. 565E+09 

310 

2 

10/19/79 

5 

47 

0 

65.1 

147.5 

115-0 

HEIGHT  OUT  OF  RANGE 

FOR  MODEL 

3 

10/19/79 

5 

47 

0 

65.1 

147.5 

116-0 

4.230E-11  9 . 366E+ 1 0 

1  .  428E+1  1 

7. 153E+1 1 

3.797E+09 

323 

2 

10/19/79 

5 

47 

0 

65.1 

147.5 

116-0 

HEIGHT  OUT  OF  RANGE 

FOR  MODEL 

3 

10/19/79 

5 

47 

0 

65. 1 

147.5 

117.0 

3.697E-11  7.953E+10 

1 . 31 OE+1 1 

6.244E+1 1 

3. 180E+09 

336 

2 

10/19/79 

5 

47 

0 

65. 1 

147.5 

117.0 

HEIGHT  OUT  OF  RANGE 

FDR  MODEL 

3 

10/19/79 

5 

47 

0 

65.1 

147.5 

118-0 

3.248E-1  1  6 . S29E+1  0 

1 . 203E+1 1 

5.477E+1 1 

2.681E+09 

349 

2 

10/19/79 

5 

47 

0 

65.1 

147.5 

118.0 

HEIGHT  OUT  OF  RANGE 

FOR  MODEL 

3 

10/19/79 

5 

47 

0 

65.1 

147.5 

119-0 

2.863E-11  5.902E+10 

1 . 106E+1 1 

4.827E+1 1 

2 . 274E+09 

362 

2 

10/19/79 

5 

47 

0 

65.1 

147.5 

119.0 

HEIGHT  OUT  OF  RANGE 

FOR  MODEL 

3 

10/19/79 

5 

47 

0 

65.  1 

147.5 

120-0 

2.545E-11  5.135E+10 

1 . 020E+ 1 1 

4.273E+1 1 

1 .9406+09 

375 

2 

10/19/79 

5 

47 

0 

65.1 

147.5 

120.0 

1.816E-11  2.71 3E+1 0 

1 . 4 1 OE+1 1 

2. 770E+1 1 

1 . 444E+09 

384 

3 

10/19/79 

5 

47 

0 

65.1 

147.5 

121-0 

2.267E-11  4.494E+10 

9.414E+10 

3.798E+1 1 

1  . 664E  +09 

388 

2 

10/19/79 

5 

47 

0 

65.1 

147.5 

121.0 

1.S89E-11  2.333E+10 

1 . 269E+ 1 1 

2.409E+1 1 

1 . 2 1 4E+09 

407 

3 

10/19/79 

5 

47 

0 

65.1 

147.5 

122.0 

2.027E-11  3.953E+10 

8.707E+10 

3.388E+1 1 

1 . 435E+09 

402 

2 

10/19/79 

5 

47 

0 

65.1 

147.5 

122-0 

1.403E-11  2 . 026E+1 0 

1 . 151 E+1 1 

2.11 5E+1 1 

1 . 03 1 E+09 

430 

3 

10/19/79 

5 

47 

0 

65.1 

147.5 

123.0 

1  .8196-11  3 . 493E+1 0 

8.068E+10 

3. 033E+1 1 

1 . 242E+09 

415 

2 

10/19/79 

5 

47 

0 

65.1 

147.5 

123.0 

1  .  249E-1  1  1.774E  +  10 

1 .051E+1  1 

1 .871E+1 1 

8 . 847E+08 

452 

3 

Figure  3.  Sample  DENMOD  Printout 
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1.1.2  Mathematical  or  Logical  Procedures 

The  main  program  first  reads  cards  1  and  2  to  determine  the  output  header, 
input  mode,  and  outputs  requested. 

If  the  time/position  input  mode  indicates  the  existence  of  an  ephemeris  file, 
card  3  is  then  read  to  give  the  earliest  and  latest  times  requested  and  the 
geophysical  parameters  to  apply  for  entire  run  (unless  APF10(1)  ^  0).  The 
dates  and  times  are  then  converted  to  modified  Julian  dates, and  fractions 
thereof,  and  stored  in  the  array  TLIM,  unless  the  months  are  0,  in  which  case 
the  values  0  and  10^  are  stored  in  TLIM.  The  first  record  on  the  file  is 
skipped  and  a  loop  is  entered.  In  each  pass  a  data  record  is  read  and  the  time 
converted  to  modified  Julian  date  and  fraction.  If  this  is  less  than  TLIM  (1), 
the  program  goes  to  the  next  record.  If  it  is  greater  than  TLIM  (2),  or  the 
end  of  the  data  is  reached,  the  program  terminates.  Otherwise  subroutine  DENS 
is  called  to  perform  the  requested  calculations,  and  the  results  are  passed  to 
subroutine  OUTDEN  for  output. 

For  card  input  the  procedure  is  similar,  except  that; 

1)  Each  card  is  processed  until  the  card  with  I  MO = - 1  is  encountered, 
causing  termination. 

2)  The  modified  Julian  date  is  not  computed  in  the  main  program,  but 
rather  in  subroutine  MODSET  via  the  call  to  DENS. 

3)  If  the  IPG  parameter  is  set,  the  main  program  sets  the  common 
variable  NLINE=0;  this  will  cause  subroutine  OUTDEN  to  start  a  new 
page. 

Subroutine  OUTDEN  is  the  output  module.  The  appropriate  data  is  passed  to  it 
through  the  common  block  /OUTDEN/.  On  the  first  call  an  initialization  phase 
is  executed.  This  includes  the  construction  of  an  execution  time  format 
statement,  which  depends  on  which  variables  are  to  be  included  in  the  printout, 
and  the  storage  of  the  total  number  of  variables  for  the  BCD  and  binary  output 
files,  if  these  latter  are  requested. 

On  each  call  OUTDEN  outputs  all  requested  variables  for  all  requested  models, 
for  a  particular  time  and  position.  A  decrementing  line  counter  NLINE  keeps 
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track  of  the  number  of  lines  left  on  the  current  page.  Setting  NLINE  to  0 
automatically  triggers  a  new  page.  This  can  be  done  by  the  main  program  since 
NLINE  is  in  the  common  block  /OUTDEN/. 


Subroutine  DENS  is  the  main  computation  subroutine,  returning  up  to  15 
variables  (mass  density,  composition,  etc.)  for  a  single  call  for  a 
particular  time,  position  and  model.  The  position  (subroutine  MODSET)  is 
converted  from  geographic  to  inertial  coordinates  using? 

2 

Rg  =  1.746647719  +  6.30038809863056d  +  0.5064  x  10-Md 

to  obtain  Rg,  the  right  ascension  of  Greenwich,  where  d  is  the  time  elapsed 
in  days  since  1  Jan.,  1950,  0  hr  UT.  Terms  due  to  nutation,  included  in  Ref. 
7,  are  neglected  here.  The  elements  of  the  Sun  are  computed  (SOL)  from  stan¬ 
dard  expressions. ^  From  these  the  solar  coordinates  are  obtained  as  follows: 

o  True  Anomoly: 

v=m  +  (2e-.25e2)sin  m  +  1.25e2  sin2m  +  1.08e2  sin3m 

where 

e  =  eccentricity 
m  =  mean  anomoly 

o  Sun-Earth  distance  (astronomical  units): 

R=(l-e2)/(l+e  cos  v) 

Mean  longitude  of  Sun  (apparent); 

L=  v  +  w  -  20.47“ /R  (0_<  L_<  2,i) 

where 

w  =  mean  longitude  of  perihelion 
o  Right  ascension: 

y  =tan_i  (tan  L  cos  E)  +  ft  •  I^(L+n/2)/nj 

where 

E  =  obliquity 

I(X)  =  largest  integer  X 

o  Declination  _ 

§  =tan_1  fsinEsin  L/^l-sin2Esin2L  1 


‘  v 

> 
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If  geophysical  values  are  not  given  directly  they  are  obtained  by  subroutine 
GETGEO  from  tables  maintained  internally  containing  up  to  100  days  of  data  read 
from  TAPE2  as  necessary  by  subroutine  APFTMP.  Instantaneous  values  required 
for  the  Jacchia  1977  model  are  obtained  by  interpolation  through  the  two 
closest  tabular  points  in  time.  The  tabular  points  for  geomagnetic  index  ap 
are  assumed  to  be  1.5  hr,  4.5  hr  etc.  on  each  day,  for  solar  activity  17  hr 
(noon  Ottawa  time).  A  similar  procedure  is  implemented  for  the  U.S.S.R.  - 
Cosmos  model  except  that  a  constant  value  period  is  assumed,  centered  at  each 
tabular  point,  during  which  the  index  is  assumed  constant  at  the  tabular  value, 
with  linear  interpolation  between  these  regions.  The  constant  value  period  is 
45  minutes  for  ap  and  6  hours  for  F10.7,  the  10.7  cm.  solar  flux.  Delay  times 
At  in  days  are  applied  to  the  time  at  which  each  index  is  computed  in  accord¬ 
ance  with  the  model  specifications: 

Jacchia  1977  U.S.S.R.  -  Cosmos 

F10.7  1.26  +  0.37sin  (H-92°)  1.0 

ap  0.1  +  0.2  cos 0.25 

where  H  is  the  solar  hour  angle  (right  ascension  at  point  -  right  ascension 
of  Sun),  and  o  is  the  geomagnetic  latitude  computed  by  the  simple  dipole 
approximation  (Ref.  1).  The  MSIS  models  do  not  require  instantaneous  values, 
simply  using  the  daily  measured  F10.7  for  the  previous  day  and  the  daily 
average  ap  (commonly  referred  to  as  Ap)  for  the  current  day. 

Smoothed  F10.7  is  assumed  constan*  for  a  given  day,  computed  as  the  average 
over  6,  6  and  3  solar  rotations  centered  on  the  given  day  for  the  U.S.S.R.  - 
Cosmos,  J77,  and  MSIS  models  respectively.  If  the  data  does  not  extend  suffi¬ 
ciently  beyond  the  day,  the  latest  available  data  is  used.  A  diagnostic  is 

printed  in  this  case.  The  average  is  corrected  for  variable  Earth-Sun  distance 

2 

R  by  multiplying  each  daily  value  by  R  prior  to  averaging,  and  dividing 

2 

the  resulting  average  by  R  . 

The  individual  models  are  then  computed  as  described  in  the  previously  mention¬ 
ed  references,  (function  RHO  for  the  U.S.S.R.  -  Cosmos  model,  subroutine  ATMDEN 
for  the  Jacchia  1977  model,  subroutine  GTS3  for  the  two  MSIS  models),  with  the 
following  notes  for  the  Jacchia  1977  model. 
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1)  The  geomagnetic  activity  thermal  component  for  each  constituent,  Alognj, 
computed  by 

4,09  "<  '  6i(^r  (°°541  *  1,  :0002(J30)2-717)  sl,’h'1(-003lT*) 

z  =  height  eAr  =  1.49 

m  =  1.7tanh  ( .005{z-100) )  =  0.00 

*N  =  0.46 

602  =  1.16  aT®  =  correction  to  exospheric  tempera- 

B0  =  0.52  ture  due  to  geomagnetic  activity 

eN£  =  1.00  Too  =  exospheric  temperature  for  zero 

6He  =  0.10  geomagnetic  activity 

2)  The  “J71  model",  Eqs  40-44  of  Ref.  1  is  used  for  the  semiannual 
variation 

3)  Atomic  Nitrogen  is  included  in  the  formulation.  Pertinent  parameters  are 

Fractional  Sea-level  Volume  =  0.00007494 
Molecular  Weight  =  14.0067 

(homopause  height  change  effect)  =  -3.70  x  10~3 
Seasonal  latitudinal  parameter  =  -0.22 

Modification  3,  above,  reflects  updates  to  the  model  by  Jacchia  since  publica 
tion  of  Ref.  1,  based  on  recent  mass  spectrometer  data^lO.  These  modifica¬ 
tions  are  discussed  in  more  detail  later. 

For  easy  reference,  the  functions  of  all  of  the  routines  shown  in  Figure  1, 
including  those  mentioned  in  the  preceding  paragraphs,  are  briefly  summarized 
in  the  following: 

DENM0D:  Input  and  control 

DENS:  Control  of  density  model  computation 

0UTDEN:  Output 

MODSET :  Conversion  of  time  and  position  to  internal  values  (modified 
Julian  date,  right  ascension,  declination) 

GETGE0:  Tabular  look  up  and  interpolation  of  solar  and  geomagnetic 
activity  indices 

APFTMP:  Construction  of  internal  tables  of  solar  and  geomagnetic  activ 
ity  indices  from  data  on  file  (TAPE2) 

Conversion  from  geomagnetic  activity  index  Kp  to  geomagnetic 
activity  index  ap 


AAP : 


AKP: 


NDTOD: 

RHO: 

READMOD : 

ATMDEM: 
GTS3: 
SOL: 
ANGLE2 : 
YRMNDY : 

NDINYR: 

DIURF: 

INTPMOD: 

BESSEL: 
FATAL : 

NEWPAGE: 


GLOBE , 

GLB: 

DCR: 

DEMS: 


Conversion  from  geomagnetic  activity  index  ap  to  geomagnetic 
activity  index  Kp 

Computation  of  the  day  number  of  the  year  from  the  date. 
U.S.S.R. -Cosmos  density  model  computation 

Initialization  for  subsequent  access  to  mass  storage  file  JSDM, 

used  for  the  Jacchia  1977  model 

Jacchia  1977  density  model  computation 

MSIS  density  model  computation 

Solar  ephemeris  computation 

Conversion  of  an  angle  to  a  value  in  the  interval  (0,  2"  ) 
Computation  of  the  calendar  date  from  the  modified  Julian 
date 

Determination  of  the  number  of  days  in  the  year  (365  or  366) 
Computation  of  Jacchia  1977  diurnal  variations. 

Tabular  look  up  of  Jacchia  1977  model  densities  and  local 
temperatures,  as  functions  of  height  and  exospheric  temperature 
Two  dimensional  interpolation 

Printout  of  fatal  error  messages  (inadequate  Jacchia  1977 
tables ) 

Skip  to  new  page  for  detailed  Jacchia  1977  printout.  This 
detailed  printout  option  is  currently  inactive,  but  may  be 
activated  by  setting  the  logical  variable  DEBUG  in  the  common 
block  GENINFO  to  .TRUE. 

Computation  of  MSIS  global,  height  independant  variations 

Computation  of  departures  from  diffusion  equilibrium  in  MSIS 
model s 

Computation  of  MSIS  height  profiles 
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1.1. 3.1  U.S.S.R.-  COSMOS 


The  U.S.S.R.-  Cosmos  model  is  an  empirical  model  for  total  mass  density  only, 
based  on  analysis  of  orbital  decay  of  Cosmos  satellites  in  the  160-300  km 
region.  The  authors  believe  the  model  to  be  valid  up  to  600  km,  as  they  have 
tried  to  relate  nighttime  and  diurnal  variations  to  the  Jacchia  19/0  model. H 
The  range  of  solar  activity  coverage  is  limited  to  values  of  75-150  in  the 
customary  units  of  10*22  w/cm^/H z.  The  principal  advantage  of  this  model 
is  its  simplicity  in  computation.  Mass  density  is  computed  directly  from 
simple  functions  without  requiring  computation  of  temperature  ana  compo¬ 
nent  densities,  as  with  the  other  models.  Extensive  tunctional  representa¬ 
tions  (as  in  the  MS15  models)  ana  tacle  look-up/  interpolation  (as  in  Jacchia 
197/  model  to  compute  diffusion  equation  solutions)  are  avoided. 

1 . 1 . 3 . 2  Jacchia  1977 

The  Jacchia  1977  ( J 7 7 )  model  is  the  first  of  the  Jacchia  series  to  include  r.a 
spectrometer  measurements  (largely  OGu-6  ana  ESRO-4).  Otherwise,  as  with 
previous  models,  the  J77  model  is  based  mainly  on  satellite  orbital  decay, 
largely  above  200km.  However  it  covers  a  broaa  range  of  solar  activity.  In 
addition  to  total  mass  density,  seven  component  densities  (0^,  0, 

N,  He,  Ar,  H),  exospheric  temperature,  and  local  temperature  are  computed. 

Geomagnetical ly  quiet  vertical  temperature  profiles  are  defined  as  functions 
of  exospheric  temperature,  which  in  turn  is  a  function  of  solar  activity, 
solar  and  geographic  declination  and  local  time.  Component  densities  are 
generated  from  this,  assuming  homogeneous  mixing  from  90km  to  lOUkm,  with 
corrections  for  0^  dissociation,  and  molecular  diffusion  above  lUOkm.  To 
account  for  observed  differences  in  local  time  phases  among  the  components, 
different  exospheric  temperatures  are  used.  Corrections  to  the  logs  of 
component  densities  are  made  for  geomagnetic  activity,  seasonal -1 ati tudi nal  , 
and  semiannual  effects.  Exospheric  and  local  temperatures  are  also  corrected 
for  geomagnetic  activity. 


A 


The  J77  model  has  the  best  coverage  over  the  various  geographical  and  geophysi¬ 
cal  conditions  and  is  particularly  suitable  for  active  conditions.  A  particu¬ 
lar  weakness  is  in  the  low  altitude  tidal  (local-time)  behavior.  It  is  known 

that  these  tidal  variations  switch  from  predominantly  diurnal  above  300km  to 

12 

predominently  semidiurnal  below  140km  ,  while  the  J77  tidal  behavior 

remains  predominantly  diurnal  at  low  altitudes.  This  is  probably  caused  by  the 
predominance  of  data  above  200km  in  the  data  base. 

1.1. 3. 3  MSI S 

The  MS  I S  models  are  based  on  in-situ  mass  spectrometer  composition  data  and 
incoherent  scatter  temperature  data.  The  satellites  range  in  height  from  160km 
to  600km.  The  data  covers  the  1966-1976  solar  cycle  (peak  solar  activity  =  150 
x  10  W/M  /Hz).  Vertical  temperature  and  composition  profiles  are 
defined  above  a  lower  boundary  at  120km  with  parameters  (exospheric  tempera¬ 
ture,  lower  boundary  temperature  and  gradient,  and  lower  boundary  composition) 
dependent  on  geographic  and  geophysical  conditions.  The  latitudi nal -tidal 
variations  are  expressed  in  spherical  harmonics  and  are  generally  considered 
an  improvement  over  the  Jacchia  representations.  The  diurnal  to  semidiurnal 
switchover  at  low  altitudes,  which  the  J 7 7  model  fails  to  represent,  is  present 
in  the  MS  I S  models.  Nevertheless  some  aspects  of  equatorial  variations  are  not 

reproduced,  possibly  a  consequence  of  heavy  weighting  of  middle  latitudes  in 
13 

the  data  base.  Geomagnetic  activity  variations  are  represented  through  a 

highly  nonlinear  function  of  the  activity  index.  However,  since  the  daily 

average  is  used,  short-term  variations  are  not  represented.  Furthermore,  J.W. 
14 

Slowey  has  pointed  out  a  near-singularity  in  this  function  at  the  geo¬ 
graphic  poles  for  zero  activity  index. 

The  MS  I S  model  with  longitude/UT  variations  (sometimes  called  MSI S  78)  is  an 
extension  of  the  original  longitude  averaged  model  by  the  addition  of  lon¬ 
gitude/UT  dependent  terms,  based  purely  on  mass-spectrometer  data,  mostly 
above  190km.  These  terms  appear  significant  mainly  in  the  polar  regions. 

Thus  this  version  may  be  more  accurate  than  the  original  MS I S  in  these  regions, 
with  no  significant  difference  at  the  equator. 


1.2 


Computational  Streamlining  of  the  Jacchia  1977  Model 


Due  to  the  increasing  complexity  of  Jacchia's  recent  density  models*  it 

becomes  necessary  to  store,  in  tabular  form,  the  solution  of  the  diffusion 

equation  for  each  constituent,  even  if  only  the  total  mass  density  is  to  be 

computed.  This  is  because  the  models  call  for  corrections  to  the  diffusion 

process  for  such  effects  as  local  time  phase  variations,  seasonal-latitudinal 

variations  and  geomagnetic  activity,  all  of  which  differ  for  the  different 

15 

constituents.  The  Jacchia  1964  model  ,  on  the  other  hand,  required 
only  temperature  corrections,  and  hence  only  the  total  mass  density  needed  to 
be  stored  as  a  function  of  height  and  exospheric  temperature.  As  a  result,  one 
must  either  set  aside  a  large  block  of  storage  (the  current  SAO  version  of 
the  Jacchia  1977  model  would  require  in  excess  of  17,000  words),  or  resort  to 
random  access  disk  storage.  The  first  approach  is  undesirable  in  almost  all 
cases,  particularly  with  small  computers.  The  second  poses  problems  in 
imp  lenient  at  ion  on  different  computers. 

Two  approaches  are  presented  here  which  take  advantage  of  the  fact  that,  with 
only  two  exceptions,  all  terms  in  the  diffusion  equation  are  analytically 
integrable.  The  exceptions  are  the  scale  height  term,  which  differs  for 
different  constituents  merely  by  a  mass  factor,  and  the  vertical  flux  term 
for  hydrogen,  which  is  important  only  for  low  altitudes.  It  therefore  is 
necessary  to  store  only  the  integral  of  the  common  mass-independent  factor  of 
the  former  term,  and  the  solution  for  hydrogen  whenever  that  constituent  is 
of  special  interest. 
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1.2.1  Condensed  Tables 


1.2. 1.1  Analysis 


The  diffusion  equation  for  the  ith  constituent,  in  the  Jacchia  1977  (J77) 
density  model  is1,16 


(1  +a. )  +  dh 
T  1  H7 


=  0 


where 

n.  =  i  constituent  number  density 
T  =  temperature 

ot j  =  i  ^  constituent  thermal  diffusion  coefficient 
H1  =  R*T/M.g 
h  =  height 


R*  =  8.31432  X  103  J (kg  -  mol)-1/deg  K 
M.=  ith  constituent  molecular  mass 
g  =  9.80665  (1  +  h /Rg ) -2m/sec2 

R  =  6.356766  X  106  meters 

6  A.  I 

=  i  constituent  vertical  flux 
D  =  diffusion  coefficient  =  2  X  10 
N  =  total  number  density 
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V  T/N~ 


The  c^.  are  assumed  values  of  -0.38  and  -0.25  for  helium  and  hydrogen  and  0 
for  all  others.  The  t>.  are  0  for  all  but  hydrogen.  Neglecting  the 

1  i  c 

vertical  flux  term  leads  to  1 

eY  E”i  F(ho-  "•  u] 

where 

F(V  h*  TJ  =  4  h  9(z)/R*T(z,  TJ  dz 
^ "  exospheric  temperature 


"i  (h-  TJ  ‘  ",  ("0.  T.) 


T(H.  TJ 
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1.2. 1.2  Tabulation 


In  practice  it  is  therefore  necessary  to  tabulate  only  F  and  the  density 
for  one  of  the  constituents  at  hQ  as  a  function  of  T«.  ,  if  h  is 
chosen  to  be  the  homopause,  100km.  If  is  chosen,  then  the  others, 
except  hydrogen,  are  given  by 

log  ni  (hQ,  T  .  )  =  log  n28  (hQ,  T  *  )  +  Q. 

where  the  subscript  i  indicates  species  by  molecular  weight  and  the  Q.  are 
constants : 

Qi  =  log  (q^./q^g)  i  1  16,  32 

Q16  =  -log  q2g  -log  ( M ’ / Mq ‘ )  +  log  [2(1-M* /M0 > )] 

Q32  =  -log  q2g  -log  (M'/M0')  +  log  M  (1  +  q  32)  *1 

V 

Where 

q.  =  sea  level  concentration  of  ith  constituent 

M*  =  mean  molecular  mass  at  100km  (Eq.  5,  ref.  1) 

M  1  =  mean  molecular  mass  at  sea  level 
o 

If  one  is  interested  in  including  the  escape  flux  term  for  hydrogen, 

special  tables  would  still  be  necessary;  however  the  storage  requirement 

would  still  be  considerably  less  than  if  separate  tables  are  used  for  all 

constituents.  Furthermore,  since  the  escape  flux  term  is  important  only 

below  500km,  it  would  be  necessary  to  store  results  only  for  that  region. 

If  one  excludes  the  escape  flux  term,  the  H  density  is  computed  using; 

h  =  500km;  eq.  17  of  ref.  1: 
o 
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log  nx  (500,  T.)  =  5.94  +  28.9  Tw1/4; 
and 

F  (500,  h,  T  „  )  =  F  (100,  h,  T  .  )  -  F  (100,  500,  T  .  ) . 

1.2. 1.3  Homogeneous  Layer  (90km  <  h  <  100km) 

For  the  homogeneous  layer  the  diffusion  equations  are  replaced  by  a 
single  barometric  equation  for  the  mass  density,  from  which  component 
densities  may  be  derived,  as  indicated  in  reference  1.  Hence  only  the 
mass  density  need  be  stored.  Alternatively  one  may  store  the  density  for 
and  derive  the  others  from  it,  as  for  the  homopause  boundary.  The 
equations  are  the  same  except  that  M'  would  be  the  mean  molecular 
weight  at  the  height  of  interest,  given  by  Eq.  5.  of  ref.  1. 

1.2. 1.4  Results 


A  subroutine  has  been  constructed  to  implement  the  above  procedures. 

Atomic  nitrogen  has  been  included,  based  on  the  most  recent  SA0 

subroutine  version,  which  incorporates  the  AE  OSS  mass  spectrometer 
9  10 

data.  *  For  this  we  have  taken  the  fractional  sea-level  composition 
given  by  Table  2,  with  the  resulting  mean-molecular  weight  = 

M  '  =  28.9586. 

Stored  in  the  subroutine  are  tables  for  the  diffusion  function  F  along  with  the 
densities  for  the  homogeneous  layer  and  the  atomic  hydrogen  solution 
including  escape  flux.  Together  these  occupy  less  storage  than  required 
for  a  single  height  segment  (4000  words)  when  the  solutions  are  tabulated 
for  each  component  on  a  random  access  file.  Hence  no  external  files  are 
required.  Run  times  on  the  CDC  6600  are  comparable  to  the  random  access 
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T 


Constituent 

Fraction  by  Volume 

Molecul ar  Mass 

Molecular  Nitrogen  (N^) 

0.78103 

28.0134 

Oxygen  (02) 

0.20953 

31.9988 

Argon  (Ar) 

0.009342 

39.948 

Hel ium  (He) 

0.000005242 

4.0026 

Atomic  Nitrogen  (N) 

0.00007502 

14.0067 

Table  2. 

Sea-Level  Composition 
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1.2.2  Analytic  Representation  of  the  Jacchia  1977  Model  Atmosphere 

Since  a  fully  analytic  solution  for  the  temperature  profiles  specified  in 
the  Jacchia  1977  (J77)  model  is  not  possible  it  seems  worthwhile  to  explore 
other  forms.  In  particular  the  Jacchia  -  Walker  (JW)profile  suggested  by 
Walker^: 


TJU<h-T„>  •  T.-<T.-V  e'°5 

leads  to  the  solution 

F(ho  h,  Tb  )  =  9eRe2 _ 

R*  (Re  +  h0)2  Ta 

where 

«  =  (h  -  h0)  (Re  +  h0)/(Re  +  h) 


T (h,  T.  ) 
Jo 


R  =  6356.766km 
e 

9e  =  gravitational  acceleration  at  Earth's  surface 
=  9.80665  m/sec2 
T0  -  T(h0.  ) 

a  =  gradient  parameter 


For  application  :o  the  Jacchia  1977  model,  it  is  prudent  to  choose  hQ  = 

125km,  since  this  is  the  inflection  point,  although  the  various  constituent 
densities  at  this  height  cannot  be  conveniently  derived  from  one  of  them,  as  at 
100  km.  The  gradient  parameter  may  be  chosen  to  fit  the  exact  J77  profile  in 
some  fashion,  such  as  to  match  the  slope  at  125km. 

Results  from  this  model  deviate  somewhat  from  the  exact  Jacchia  1977 
model.  In  particular,  for  an  exospheric  temperature  of  900K  the  local 
temperature  deviation  is  20K  at  160km  and  the  mass  and  argon  densities 
differ  by  6%  and  12%  at  225km. 

To  obtain  better  fits  the  form 
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has  been  explored,  where  the  function  f  is  chosen  to  retain  the  closed-form 
integrability  of  the  diffusion  equation.  In  addition  it  is  desirable  that  the 
resulting  total  temperature  profile  fit  smoothly  the  exact  J77  profile  at 

h  =  h0  =  125km  and  maintain  the  inflection  point  character  (d^T/dh^  =  0). 

It  is  evidently  the  lack  of  this  latter  characteristic  which  makes  the  simple 
JW  form  undesirable  and  possibly  accounts  for  much  of  the  error  resulting  from 
its  use.  These  modifications  lead  to  the  following  boundary  conditions  at 

S=  o- 

f  =  df/dC  =  0 

2  2  2  2 

d  T/d-  =  -G  r/T  c 
o  o 

where 


V^h  -  ho 

+  2/(Re  +  hQ)  =  G0/(Tot-  Tq)  +  2/(Re  ♦  h0) 


For  convenience  the  functional  form  Is  chosen  to  have  an  exponentially 
decreasing  character  with  parameters  hopefully  related  too  .  One  such 
form  is 


f  = 


C/(  2S2Tn 

i  0 


f8  (5) 


where 


and  8  is  an  adjustable  parameter.  To  provide  further  flexibility  the 
form 


f(f)  =  J 


V> 


fg  U)  ♦  c/,fs  (c) 

l  2 


has  been  chosen.  The  additional  term  j  =  C  £f Q  satisfies 

2  P 

2 


J  “  j  ■  j' 


■S 
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at  C  *  0. 


For  this  profile  the  solution  is  given  by 


p(h0,  h,  TJ  =  FJW  +  a.  q|  ♦  l/3)/B 

+C  I.L  (B  £a .  -  b.)  q  1  + 

2  P~ 1  2  I  1  2 


where  FJW 


=  Jacchia-Walker  solution  (for  TJW) 
Cj  =  -Gq  5/(2 
a^  -  ”^2  ~  - 1 >  “1/2 


bi  =  a  i  /  i 

K  =  geRe2/[R*(Re  +  hQ)2]. 


i 

li/iaU2} 
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1.2.2. i  Results 


18  19 

The  Fletcher-Powel 1  non-linear  function  minimization  method  ’  has 
been  used  to  determine  the  parameters  8,6,  c  t0  minimize  the  sum  of 
the  squared  residuals  1  2  2 

Q  ■  I(T  -  Tj 77) 2 


for  T“=  500,  700,  900,  1100,  1300,  1500,  1700,  1900.  By  inspection  the 
results  are  found  to  be  reasonably  well  represented  by: 


6 

1 

6 

2 
C 

2 


1 


2o  1  +  10 


T„  S  1100K 

(T.  - 


0.0215  -  0.005  (T^  -  500)/ 200 

10~s(0.0566  -  0.08X  +  0.04X2);  X  =  (T^  -  900)/200 


1100K  <  T  <  1500K 

OO 

B  =  0.0385  -  0.012y  +  0.0123y2;  y  *  (T^-  1100)/400 

6z  =  0.0065  +  0 ,0167y 

C2  =  10's(0.0166  -  0.4548y2) 

T  >  HOOK 
00 

B  =  0.0388  -  0.0045z;  z  =  (T^  -  1500)/400 

6  =  0.0232  -  0.0040Z 
2 

C  =  -10'5  (0.4382  +  0.0387z) 

2 
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Mass  and  constituent  densities  were  computed  using  the  Jacchia-Bass  (JB) 
formulation  presented  here  for  the  T,*  region  500K  -1900K  at  1Q0K  steps  and  the 
height  region  130km  -  1000km  at  10km  steps.  The  results  were  compared  with 
those  from  the  exact  J77  temperature  profiles,  except  that  vertical  flux  is 
ignored  for  H.  Table  3  summarizes  maximum  %  deviations  for  a  6-constituent  gas 
(02,  0,  N2,  He,  Ar,  H)  with  the  total  mass  density  given  by 


o 


Z 

i  =  l 


Mi  nt/A; 


A  =  Avoqadro's  Number 


It  should  be  noted  that  a  true  J77  calculation  would  call  for  different  exo¬ 
spheric  temperatures  for  the  different  constituents  to  model  the  different 
diurnal  phases;  here  all  the  n^  are  computed  for  the  same  T».  Argon, 
with  the  largest  mass,  yields  the  largest  disagreement  among  the  con¬ 
stituents.  The  agreement  is  quite  good,  well  within  typical  model -experiment 

20 

differences  such  as  those  published  by  Forbes,  Marcos  and  Gillette  ,  and 
21 

Sharp  and  Prag.  It  should  be  noted  that  the  maximum  deviations  are  all 
within  the  1100K  -  1500K  region,  which  proved  to  be  the  most  difficult  to 
fit.  Outside  this  region  all  JB  mass  densities  agree  with  J77  within  1%  at 
all  heights.  Maximum  temperature  deviations  are  13K  at  T„  =  1300K,  h  = 

350km.  Outside  the  1100K-1500K  region  the  maximum  disagreement  is  5K. 


1.2. 2. 2  Conclusions 


It  has  been  shown  that  analytically  solvable  temperature  profiles  can  be 
adjusted  to  realistically  represent  temperature  -  height  variations  above 
12bKm,  allowiny  one  to  reduce  tabular  storage  requirements  and/or  avoid  the 
inconvenience  of  interpolation.  It  is  therefore  recommended  that  similar 
profiles  be  used  in  future  atmospheric  density  modelling.  The  possibility  of 
developing  analytically  integrable  profiles  below  125Km  will  be  studied  at  a 
later  date. 
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'  1 


TOTAL  DENSITY  Ar  DENSITY 


h 

300Km 

1.1 

(h  =  200,  T  0=  =  1200) 

1.5 

(h  =  190,  T.  = 

1200) 

h 

500Km 

2.1 

(h  =  500,  T»  =  1300) 

5.6 

(h  =  500,  T-  « 

1300) 

h 

lOOOKm 

3.0 

(h  =  750,  T.  =  1300) 

8.5 

(h  =  1000,  Too 

=  1300) 

Table  3.  Maximum  Absolute  Value  of  JB-J77  Deviations  (%) 


37 


1.3.  Satellite  Orbit  Prediction  Error  Evaluations 


Density  Model  evaluatior  by  means  of  satellite  orbit  prediction  errors  has 

continued  as  part  of  an  on-going  effort  to  assist  tracking  and  surveillance 

agencies  in  selection  of  models  for  performance  of  these  functions.  Although 

models  also  are  evaluated  in  terms  of  their  inherent  accuracy  compared  with 
20  21 

existing  data  ’  ,  prediction  accuracy  evaluation  should  provide  a  more 

direct  means  of  assessing  a  model's  use  in  tracking.  Calibration  could  be 
established  between  prediction  accuracy  and  inherent  accuracy.  Furthermore 
prediction  accuracy  evaluation  can  be  used  to  emphasize  those  variables 
important  to  a  particular  satellite  of  interest,  such  as  geomagnetic  acti¬ 
vity,  while  de-emphasizing  less  important  variables  such  as  local  time  for  a 
polar-orbiting  satellite. 

The  principle  software  vehicles  for  these  studies  are  Doppler  beacon  analysis 
program  CELEST^  and  the  radar  data  analysis  package  CADNIP/BADMEP^. 

This  report  will  focus  primarily  on  CELEST,  although  it  substantially  applies 
to  CADNIP/BADMEP  as  wel 1 . 

1.3.1  Procedures 


The  basic  procedures  have  been  described  in  previous  reports  (for  example. 

Ref.  6).  In  summary,  two  CELEST  runs  are  required  per  case.  In  the  first  run, 
data  is  fit  over  a  selected  span  and  the  resulting  trajectory  extended  by 
numerical  integration  over  a  subsequent  predict  span.  In  the  second  run 
data  is  fit  over  a  span  overlapping  the  first  fit  span  and  including  data  in 
the  predict  span  in  order  to  obtain  a  "true"  trajectory  to  compare  with  the 
predicted  trajectory  obtained  from  the  first  run.  In  each  run  the  fit  is 
obtained  by  a  least  squares  adjustment  of  Keplerian  elements  at  the  start  of 
the  fit  span,  and  the  drag  coefficient  for  the  whole  fit  span,  the  latter 
of  which  essentially  acts  as  a  scale  factor  to  account  for  average  error  in  the 
chosen  density  model.  This  procedure  is  repeated  for  several  density  models  and 
time  periods.  Statistics  (mean  magnitude  error,  mean  algebraic  error,  standard 
deviation)  are  then  taken.  Automated  procedures^  are  used  to  simplify 
operations. 


3B 


Time  periods  are  chosen  to  reflect  as  nearly  as  possible  overall  variations 
in  geomagnetic  activity  for  low  altitude  satellites.  Fit  spans  of  one  day 
followed  by  15  hour  prediction  periods  are  generally  sufficient  to  establish 
trends.  Longer  time  spans  are  needed  for  higher  altitude  satellites.  Care 
is  taken  to  avoid  orbit  adjusts. 

1.3.2  Program  Changes 

The  MSIS  78  model  has  been  added  to  CADNIP/BADMLP ,  replacing  the  Harris- 
Priester  density  model.  Since,  as  previously  noted,  this  model  differs  from 
the  MSIS  model  only  in  terms  based  mainly  on  data  above  190Km,  the  model  has 
not  been  added  to  CELEST.  Table  4  indicates  the  various  models  available  in 
these  programs.  The  NWL  model  is  described  in  Ref.  24;  and  the  Forbes-Garrett- 
Gillette  model,  in  Ref.  25.  Other  models  not  previously  discussed  are  described 
in  Ref.  23  and  references  listed  therein. 

Calculation  of  solar  and  geomagnetic  activity  indices  has  been  upgraded  to 
include  the  linear  interpolation  scheme  described  in  section  1.1.2  for  the 
Jacchia  1977  model  in  program  DENMOD.  Previously  these  functions  had  been 
regarded  as  constant  for  the  time  segment  of  each  tabular  value  (1  day  for 
solar  activity,  3  hrs  for  geomagnetic  activity),  jumping  di scontinuously  at 
the  beginning  of  each  new  time  segment.  The  linear  interpolation  method  more 
accurately  reflects  the  processing  used  in  the  development  of  the  Jacchia 
model s . 

Calculation  of  smoothed  solar  activity  has  been  modified  to  decouple  the 
variation  in  solar  distance  over  the  smoothing  period,  also  as  described  in 
Section  1.1.2.  The  reason  for  this  procedure  is  that  the  smoothed  solar 
activity  is  supposed  to  represent  the  solar  disk  component  (in  contrast  to  the 
active-region  component)  of  the  solar  radiation  striking  the  Earth's  surface  on 
the  day  for  which  the  smoothing  is  taken,  vis.,  the  mid-point  of  the  smoothing 
period.  Underlying  this  approach  is  the  fact  that  the  average  (smoothing)  of 
the  daily  index  at  a  fixed  point  in  space  over  a  period  of  time  is  known  to 
correlate  with  the  disk  component  of  the  solar  radiation  reaching  that  point  at 
the  middle  of  this  time  period.  The  daily  activity  index  is  the  10.7cm  flux 
measured  at  the  Earth's  surface,  and  therefore  it  varies  inversely  as  the 
square  of  the  the  Sun-Earth  distance  R.  Therefore  each  value  must  be 
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Model  Number 

CELEST 

CADNIP/BADMEP 

0 

NWL 

1 

Jacchia  1964 

Jacchia  1964 

2 

1966  Supplements 

1966  Supplements 

3 

Jacchia  1971 

Jacchia  1971 

4 

U.S.S.R.  -  Cosmos 

U.S.S.R.  -  Cosmos 

5 

Jacchi a-Walker-Bruce 

Jacchi a- Wal ker-Bruce 

6 

Jacchia  1977 

Jacchia  1977 

7 

Lockheed/NASA 

Lockheed/NASA 

8 

MSI  S 

MSI  S 

9 

1962  U.S.  Standard 

1962  U.S.  Standard 

10 

MSI S  78 

11 

DENSEL 

12 

Jacchia  1970 

Jacchia  1970 

13 

Jacchia  1973 

Jacchia  1973 

14 

Forbes-Garrett- 

Gillete  Model  B 

Table  4.  Density  Models  Available 
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multiplied  by  the  daily  r2  value  to  obtain  the  flux  at  a  fixed  distance, 

1  A.U.  The  resulting  average  is  then  divided  by  R^  for  the  day  in  question 
to  give  the  correct  value  at  the  Earth's  surface  for  that  day.  If,  for  example, 
for  the  day  in  question  the  Earth  is  at  aphelion,  failure  to  implement  this 
procedure  would  result  in  an  overestimate  of  the  disk  component  of  the 
total  solar  flux  striking  the  Earth's  surface  that  day. 

1.3.3  Results  and  Interpretation 

Of  the  three  statistical  measures  mentioned  previously,  mean  absolute  value 
error  and  standard  deviation  are  generally  accepted  as  the  best  criteria  for 
comparison.  The  algebraic  mean  is  used  only  to  detect  possible  bias  in  the 
model  for  the  selected  sample.  The  results  are  partitioned  into  high  and  low 
geomagnetic  activity  bins  to  permit  determination  of  which  models  are  best 
under  particular  conditions  (other  variables  can  be  partitioned  if  appro¬ 
priate).  Overall  averages  are  then  corrected  for  incorrect  weighting  with 
respect  to  geomagnetic  activity  or  other  variables. 

I .4  Atmospheric  Density  Determination 

In-situ  measurement  by  accelerometers,  mass  spectrometers,  and  ion  density 
gauges  has  replaced  satellite  orbit  decay  analysis  as  the  primary  means  of 
measuring  atmospheric  density.  These  instruments  provide  time  and  spatial 
resolution  unavailable  from  orbital  decay  analysis  because  of  the  inherent 
smoothing  present  in  the  latter.  In  addition  mass  spectrometers  and  ion  gauges 
provide  composition  data.  Nevertheless  orbital  decay  methods  should  be  useful 
for  the  following  reasons: 

1)  Supplement  to  in-situ  data  base 

2)  Calibration  for  in-situ  measurements  by  comparision  of  orbit- 
averaged  or  orbit-integrated  densities 

3)  Diagnostic  tool  when  in-situ  data  unavailable,  for  instance  to 
determine  relation  of  prediction  errors  to  density  variation 


Density  determinations  have  thus  continued,  using  CELEST  as  described  in  Ref. 
24  in  conjunction  with  Doppler  beacon  data  received  from  the  Defense  Mapping 
Agency  (DMA). 

1.4.1  Data  File  Conversion 


Critical  to  the  operation  of  CELEST  is  the  proper  creation  of  certain  data 
files  (TAPE14:  preprocessed  observations;  TAPE 19 :  master  Sun-Moon/coordinate 
transformation  tables)  generated  elsewhere.  This  data  is  sent  from  DMA  in 
BCD  format  and  must  be  converted  to  binary  for  use  by  CELEST.  Formats  for 
this  data  are  given  in  Tables  5-8. 
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RECORD 

TYPE 

WORD 

NUMBER 

VARIABLE 

TYPE 

VARIABLE 

NAME 

FORMAT  USED 
FOR  ENCODING 

1 

1 

A 

WORD 

A6 

2 

I 

SAT 

14 

2 

1 

R 

YEAR 

F6.0 

2 

R 

DAY 

F6.0 

3 

D 

TTCA 

E20.14 

4 

A 

ISTA 

A6 

5 

I 

ICLAS 

12 

6 

I 

ITYPE 

12 

7 

D 

OCTOL 

E20.14 

8 

I 

NO 

13 

9 

I 

NBI 

13 

10 

R 

TEMPR 

F5.0 

11 

R 

PRESS 

F5.0 

12 

R 

HUM 

F5.0 

13 

I 

ITI 

12 

14 

I 

IF 

12 

15 

I 

IP 

12 

16 

I 

IQPR 

16 

3 

1 

D 

XL  AM 

E20.14 

2 

D 

PHI 

E20.14 

3 

D 

ALT 

E20.14 

4 

D 

FS 

E20.14 

5 

R 

DUMMY  (1) 

F10.0 

6 

R 

DUMMY  (2) 

F10.0 

7 

R 

DUMMY  (3) 

F10.0 

8 

R 

DUMMY  (4) 

F10.0 

4 

1 

D 

DA ( 1 , 1 ) 

E19.13 

2 

D 

DA (1,2) 

E 19 .13 

3 

D 

DA ( I ,3) 

E12.6 

4 

D 

DA(I ,4) 

E10.4 

5 

D 

DA (1+1,1) 

E19.13 

6 

D 

DA ( I +1,2) 

E19.13 

7 

D 

DA( I+l ,3) 

EI2.6 

8 

D 

DA (1+1,4) 

E10.4 

Record  type  4  is 

repeated  until 

all  data  points 

are  written 

(1=1, NO) 

Record  types  2,  3 

and  4  are  repeated  for  each  data  pass. 

A  description  of 

each  of  the  above  variables  is 

contained  in 

Table  6. 

Table  5.  CELEST 

Time  Corrected 

Observations : 

BCD  Data  File  Description 
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Record 

Word 

Word 

Type 

Number 

Type 

Quantity 

1 

1 

A 

10BC 

2 

I 

SAT 

2 

1 

R 

YEAR 

2 

R 

DAY 

3 

R 

TTCA 

4 

A 

ISTA 

5 

I 

ICLAS 

6 

1 

ITYPF 

7 

R 

OCTOL 

8 

1 

NO 

9 

I 

NBI 

10 

R 

TEMPR 

11 

R 

PRESS 

12 

R 

HUM 

13 

I 

ITI 

14 

I 

IF 

15 

I 

IP 

lb 

I 

IQPR 

17 

R 

XL  AM 

18 

R 

PHI 

19 

R 

ALT 

20 

R 

FS 

21-24 

R 

DUMMY ( 1 ) -DUMMY (4 ) 

3 

1 

R 

DA(1 ,1 ) 

2 

R 

DA(1,2) 

3 

R 

DA( 1 ,3 ) 

4 

• 

R 

DA( 1 ,4 ) 

• 

• 

4*N0 

R 

• 

DA(N0,4) 

1  type  1  record  per  file 
1  type  2  record  and  1  type  3 
record  per  data  pass 


Table  6.  CELEST  Time  Corrected  Observations:  Binary 


Description 
"QBS  FI" 

NWL  Satel 1 ite  No. 

Observation  year 
Observation  day 
Predicted  time  of 
closest  approach 
(se-  from  midnight) 
Station  number 
Data  Class 
Data  type 
0-C  tolerance  for 
filtering  data 
Number  of  obs. 

Unused  (=0) 
Temperature 
Pressure 
Humidity 
Unused  (=0) 

Unused  (=0) 

Unused  (=0) 

Q- number 

Station  longitude 
(Deg) 

Station  latitude 
(Deg) 

Station  altitude 
(km) 

Satellite  frequency 
Unused  (=0) 

Time  of  1st  obs. 
Observation  value 
Sigma  for  obs. 

Tag  (0=good  obs) 


Data  File 
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Record  type 

Format 

#Words 

1 

3  E20.8 

3 

2 

5  E18.8 

20 

3 

5  E18.8 

20 

Record  type  1  occurs  once  per  year 

Record  type  2  ocurs  once  per  day,  4  lines  each 

Record  type  3  occurs  once  per  year,  following  daily  Record  type  2 
data  for  that  year 

Data  described  in  Table  8 


Table  7.  CELEST  Master  Sun-Moon/Coordinate  Transformations 

BCD  File  Description 


Record 

Word 

Word 

Type 

Number 

Type 

Quantity 

1 

1 

R 

RTYPE 

2 

R 

RYR 

3 

R 

ROD 

2 

1 

R 

RTYPE 

2 

R 

DRDA 

3 

R 

XSUN 

4 

R 

YSUN 

5 

r> 

ZSUN 

6 

R 

XM00N  1 

7 

R 

YM00N  1 

8 

R 

ZM00N  1 

9 

R 

XM00N  2 

10 

R 

YM00N  2 

11 

R 

ZM00N  2 

12 

R 

DLSI 

13 

R 

EP 

14 

R 

DELEP 

15 

R 

DLT 

16 

R 

0LH 

17 

R 

P 

18 

R 

Q 

19,20 

1 

R 

RTYPE 

2-20 


Description 

Any  positive  number 
Year  (last  2  digits) 

Number  of  days  in 
the  year 

Any  positive  number 
Day  number  (1  Jan  =1.0) 
Rectangular  coordinates 
(km)  of  sun  in  1950 
inertial  reference 
system  at  midnight 
(GMT)  of  day 

Rectangular  coordinates 
(km)  of  moon  in  1950 
inertial  ref.  system 
at  midnight  (GMT)  of  day 

Same  as  words  6-8, 
except  at  noon 

Nutation  in  longitude 
(radians) 

Mean  obliquity  (radians) 
Nutation  in  obliquity 
(radians) 

Seasonal  correction  to 
Earth's  rotation  (sec.) 
Equation  of  equinoxes 
(radi ans ) 

Polar  motion  correction 
parameters-added  to  2nd 
and  1st  parameters, 
respectively  on  RCC08  card; 
usually  =  0 
Unused  (=0) 

Negative  number,  indicates 
this  is  last  record  of  year 

Unused 


Words  12-18,  record  type  2  are  at  midnight,  GMT. 


Table  8.  CELEST  Master  Sun-Moon/Coordinate  Transformation:  Binary  File  Description 
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Programs  DMABIN  and  BINAR  perform  the  required  conversion  of  the  observation 
and  Sun-Moon/transformation  data,  respectively,  from  the  BCD  to  binary  formats 
required  by  CELEST.  The  input  tape  in  each  case  has  the  following  specifi¬ 
cations: 

Logical  file  name:  TAPE  1 
Density:  HI  (556  BPI,  7-track) 

Record  Manager:  RT=S,  BT=C 

The  output  tape  in  each  case  is  TAPE  2.  Program  DMABIN  requires  the  follow¬ 
ing  input  card: 


Col 

Variable 

Format 

Description 

1-5 

NFILE 

15 

Number  of  files 

on  input  tape 

6-10 

IREW 

15 

l=rewind  at  end  of 

run 

0=don‘t  rewind 

The  observation  input  tapes  generally  contain  1-5  files  (days)  of  data  each, 
but  conversion  to  binary  and  selection  of  higher  density  permit  much  more  to  be 
written  on  one  output  tape.  Thus  several  runs  of  DBAMIN  may  be  executed  in 
one  job,  replacing  input  tapes  between  runs.  Thus  IRE W  should  be  0  for  all 
but  the  last  run,  when  IREW  should  be  1  to  prevent  the  system  from  attempting 
to  read  blank  tape  in  subsequent  use.  The  output  tape  would  thus  have  an 
end-of-file  after  each  file  (day)  of  data  and  be  terminated  by  proper  end-of- 
infonnation  markers  if  rewound  at  end  of  last  run. 

Program  BINAR  requires  no  input  other  than  the  tape  which  usually  contains 
two  years  of  data.  In  practice  the  output  can  be  catalogued  as  a  permanent 
file  on  disk;  thus  rewind  is  not  necessary. 

1.4.2  Operational  Procedures 


Operation  of  CELEST  is  as  described  in  Ref.  24  except  that,  as  noted  previously, 
observations  received  from  DMA  are  already  preprocessed.  Thus  execution  of  the 


preprocessor  module  in  CELEST  is  skipped.  The  data  span  (usually  one  day  and 
the  last  4  1/2  hrs  of  the  previous  day  for  overlap)  can  be  divided  into  several 
drag  segments.  The  Keplerian  elements  at  the  start  of  the  fit  span  plus  a  drag 
coefficient  for  each  drag  segment  are  adjusted  to  perform  a  least-squares  fit 
to  the  data.  Density  at  perigee  can  be  estimated  as 

D  =  (<VCTH)  ^ m0del 

where 

C.j  =  drag  coefficient  for  the  segment  in  which  the  perigee  occurs 
CjH  =  tf3oretical  drag  coefficient  (=2.2) 

Dmodel  =  Densit-y  at  perigee  for  the  selected  model 

1.4.3  Rapid  Density  Variations 

Very  short  drag  segment  durations  are  desirably  to  obtain  the  time  resolution 
required  to  study  the  rapid  variations  which  occur  during  magnetic  storms. 
However,  tracking  data  distribution  in  time  and  geopotential  modelling  errors 
place  limits  on  the  time  resolution  which  can  be  obtained.  The  distribution 
and  number  of  station  passes  in  one  revolution  are  highly  variable,  introducing 
much  noise  to  the  results  if  such  short  segments  are  used.  Geopotential 
modelling  errors  begi-  to  be  absorbed  significantly  into  the  drag  coefficient 
determined  by  least  squares  for  segments  smaller  than  four  revs,  although  this 
effect  may  still  be  dominated  by  real  drag  variations  if  the  latter  are  suf¬ 
ficiently  large.  Thus  two  revolutions  appears  to  be  the  best  possible  time 
resolution  that  can  be  attained.  Density  variations,  however,  may  occur  on  a 
much  smaller  time  scale  during  high  geomagnetic  activity,  because  of  the  very 
sharp  latitudinal  dependence  which  may  occur.  Therefore  perigee  density,  as 
defined  above,  is  not  a  reliable  or  meaningful  quantity.  The  most  meaningful 
result  is,  rather,  a  suitable  average  of  drag  or  density  over  the  duration  of 
the  drag  segment. 

For  this  reason,  CELEST  has  been  augmented  with  the  capability  to  integrate 
numerically  the  drag  acceleration  magnitude  over  a  segment  to  obtain  the 
average  drag  over  a  segment  by: 

=  cd  y*1,  9d  (t)  dt/^'  -  *») 
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where 


Cq  =  drag  coefficient  determined  for  the  segment 
t0  =  start  time  of  segment 
t,  =  end  time  of  segment 
9o(t)  =  ADV2/2M 

A  -  satellite  area  to  mass  ratio 
M 

D  =  atmospheric  density  according  to  chosen  model 
V  =  satellite  velocity 

The  function  go(t)  is  integrated  using  the  Gauss  -  Jackson  formulation 
discussed  in  Reference  22,  expanding  the  3  dimensional  position  vector  x 
being  solved  for  to  a  4  dimensional  vector: 

X  =  G  (x,  x,  t) 

where  the  first  three  components  of  x  are  the  original  Cartesian  position 
components  and  the  fourth  component  corresponds  to  the  drag: 

64  =  9D 

Thus  X4  is  the  result  of  interest.  Initial  conditions  are 


X4  =  X4  =  0. 


The  average  acceleration  over  a  drag  segment  is  then: 
Oq  =  [X4(t. )  -  X4  (to)]/(t,  -  to) 


The  orbit  generation  module  of  CELEST  was  upgraded  to  perform  this  expanded 
computation. 

The  average  drag  predicted  by  the  chosen  model  is  computed  for  comparison 
(Figures  4  and  5)  simply  by  setting  Cq  to  a  fixed  value.  These  results 
clearly  indicate  a  greater  response  of  the  model  to  geomagnetic  activity 
increases  for  the  particular  active  period  shown  than  the  actual  atmospheric 
response  as  deduced  from  the  tracking  data. 


49 


ar<]  4- 

or 


CJlSI 
<_i  ii 

aud 
Dl- 
OwCt 
a  a 
<fc_i 
aujoa 
oo 

LUO 

(31 CO 
cr  *- 


OiUJt- 
00(0 
O  UJ 


r-uj 

UJ 

a: 

X 


to 


(  S3330b3 I W )  N0Ilbd3 1333b  iyH9-l)d« 


Figure  4.  Comparison  of  3  Hour  Averaged  Drag  Acceleration 
predicted  by  J77  Density  Model  (Triangles)  with 
that  deduced  frc  Satellite  Trackin'  Data  (Crus 


1.5 


References 


1.  Jacchia,  L.G.,  "Thermospheric  Temperature,  Density,  and  Composition:  New 
Models",  Smithsonian  Astrophysical  Observatory,  Special  Report  Number  375, 
1977. 

2.  Hedin,  A.E.,  Salah,  J.E.,  Evans,  J.V.,  Reber,  C.A.,  Newton,  G.P., 

Spencer,  N.W.,  Kayser,  D.C.,  Alcayde,  D.,  Bauer,  P.,  Cogger,  L.,  and  McClure, 
J.  P.,  "A  Global  Thermospheric  Model  Based  on  Mass  Spectrometer  and  Incoherent 
Scatter  Data-MSIS  1.  N^  Density  and  Temperature",  J.  Geophys  Res.  82,  p. 

2139,  1977;  Hedin,  A.E.,  Reber,  C.A.,  Newton,  G.P.,  Spencer,  N.W.,  Brinton, 

H. C.,  Mayr,  H.G.,  and  Potter,  W.E.,  "A  Global  Thermospheric  Model  Based 

on  Mass  Spectrometer  and  Incoherent  Scatter  Data  -  MSIS  2.  Composition",  J. 
Geophys.  Res.  82,  p.  2148,  1977. 

3.  Kayser,  D.C.,  Breig,  E.L.,  Power,  R.A.,  Hanson,  W.3.  and  Nier,  A.O., 
"Direct  In  Situ  Measurements  of  Thermospheric  Temperature",  J.  Geophys.  Res. 
84A,  p.  4321,  1979. 

4.  Elyasberg,  P.E.,  Kugaenko,  B.  V.,  Synitsyn,  V.  M.,  and  Voiskovsky,  M.I., 
"Upper  Atmosphere  Density  Determination  from  Cosmos  Satellite  Deceleration 
Results",  Space  Research  XII,  p.  727,  1972. 

b.  Hedin,  A.  E.,  Reber,  C.  A.,  Spencer,  N.W.,  Brinton,  H.C.,  and  Kayser, 

D.C.,  "Global  Model  of  Longitude/UT  Variations  in  Themospheric  Composition 
and  Temperature  Based  on  Mass  Spectrometer  Data",  J.  Geophys.  Res.  84A,  p. 

I. ,  1979. 

6.  Logicon,  Inc.,  "Research,  Analysis,  Development,  and  Application  of 
Highly  Integrated  Systems  of  Multi-Phases  of  the  Physics  of  the  Upper  Atmos¬ 
phere",  Final  Report,  AFGL-TR-78-0204,  1978. 

7.  Minka,  K.,  "Orbit  Determination  and  Analysis  by  the  Minimum  Variance 
Method",  AFGL-65-b79 ,  p.  12,  1965. 

8.  "Explanatory  Supplement  to  the  Astronomical  Ephemeris  and  the  American 

Ephemeris  and  Nautical  Almanac",  Her  Majesty's  Stationery  office,  publ., 
London,  p.  98,  1961. 


52 


9.  Mauersberger ,  K.,  Engebretson,  M.J.,  Potter,  W.E.,  Kayser,  D.C.,  and 
Nier,  A.O.,  "Atomic  Nitrogen  Measurements  in  the  Upper  Atmosphere",  Geophys. 
Res.  Let.  2,  p.  337,  1975. 

10.  Engebretson,  M.J.,  Mauersberger,  K.,  and  Potter,  W.E.,  "Extension  of 
Atomic  Nigrogen  Measurement  into  the  Lower  Thermosphere",  J.  Geophys.  Res.  82, 
p.  3291,  1977. 

11.  Jacchia,  L.G.,  "New  Static  Models  of  the  Thermosphere  and  Exosphere 
with  Empirical  Temperature  Profiles"  Smithsonian  Astrophysical  Observatory, 
Special  Report  Number  313,  1970. 

12.  Sharp,  L.R.,  Hickman,  D.R.,  Rice,  C.J.,  and  Strauss,  J.M.,  "The  Altitude 
Dependence  of  the  Local  Time  Variation  of  Thermospheric  Density"  Geophys.  Res. 
Let.  j>,  p.  261  ,  1978. 

13.  Forbes,  J.M.,  and  Marcos,  F.A.,  "Tidal  Variations  in  Total  Mass  Density 
as  Derived  from  the  AE-E  Mesa  Experiment",  J.  Geophys.  Res.  84A,  p.  31,  1979. 

14.  Slowey,  J.W.,  private  communication. 

15.  Jacchia,  L.G.,  "Static  Diffusion  Models  of  the  Upper  Atmosphere  with 
Empirical  Temperature  Profiles",  Smithsonian  Astrophysical  Observatory, 

Special  Report  Number  170,  1964. 

16.  COSPAR  Working  Group  IV,  "COSPAR  International  Reference  Atmosphere  1965 
(CIRA  1965)",  Noth-Holland  Publishing  Company-Amsterdam,  P.  7,  1965. 

17.  Walker,  J.C.G.,  "Analytic  Representation  of  Upper  Atmospheric  Densities 
Based  on  Jacchia's  Static  Diffusion  Models",  J.  Atm.  Sci .  22,  p.  462,  1965. 

18.  Fletcher,  R.,  and  Powell,  M.J.D.,  "A  Rapidly  Convergent  Descent  Method 
for  Minimization",  Computer  Journal  6^,  p.  163,  1963. 

19.  International  Business  Machines  Corp. ,  "System/360  Scientific  Sub¬ 
routines  Package",  Fifth  Edition,  p.  221,  1970. 


* 


53 


20.  Forbes.  J.M.,  Marcos,  F.A.,  and  Gillette,  D.F.,  "An  Evaluation  of 
Thermospheric  Models",  AFGL-TR-78-0140 ,  1979. 

21.  Sharp,  L.R.,  and  Prag,  A.B.,  "The  Construction  of  Thermospheric  Density 
Correction  Tables  for  Use  in  Satellite  Trajectory  and  Lifetime  Predictions", 
SAMS0-TR-78-87 ,  1978. 

22.  O'Toole,  J.W.,  “Naval  Surface  Weapons  Center  Reduction  and  Analysis  of 
Doppler  Beacon  Satellite  Receivers  Using  the  CELEST  Computer  Program", 
Satellite  Doppler  Positioning,  Volume  2  of  Proceedings  of  the  International 
Geodetic  Symposium,  October,  197b. 

23.  Bramson,  A.S.,  and  Slowey,  J.W.,  "Some  Recent  Innovations  in  Atmospheric 
Density  Programs",  AFCRL-TR-74-0370,  1974;  "Atmospheric  Model  Evauation", 
AFCRL-71-0543 ,  1971  . 

24.  Bass,  J.N.,  Bhavnani  ,  K.H.,  and  Hussey,  I. A.,  "Atmospheric  Density 
Determination  from  Analysis  of  Doppler  Beacon  Satellite  Data",  AFCRL-TR-75- 
0176,  1975. 

25.  Garrett,  H.B.,  "An  Updated  Empirical  Density  Model  for  Predicting  Low- 
Altitude  Satellite  Epheinerides ,"  AFCRL-TR-75-0158 ,  1975. 


54 


2.0  Satellite  Ephemerides 

Critical  to  the  use  of  satellites  for  geophysical  measurement  and  exploration 
tasks  is  the  need  to  have  knowledge  of  vehicle  position  as  a  function  of  time. 
In  some  appl ications,  the  requirements  for  accuracy  can  be  very  exacting.  In 
general  terms,  the  problem  involves  some  combination  of  smoothing,  inter¬ 
polating  and  extrapolating  radar  observations.  This  section  is  devoted  to 
discussion  of  LOKANGL,  a  program  that  has  proven  of  great  usefulness  in 
furnishing  orbital  information  to  support  a  variety  of  satellite  projects.  A 
system  of  coordinates,  known  as  solar/magnetospheric,  has  been  found  useful 
in  the  study  of  solar  interactions  with  the  magnetosphere.  In  support  of  the 
SCATHA  effort,  the  capability  has  been  added  to  Program  LOKANGL  to  express 
satellite  ephemerides  in  this  system  and  to  plot  the  results. 

2.1  LOKANGL 


2.1.1  Introduction 

Program  LOKANGL  is  a  utility  program  for  general  calculations  of  satellite 
ephemerides.  It  has  assumed  this  role  with  the  general  availability  of  high 
quality  orbital  elements,  which  form  the  major  input  to  the  calculations. 

The  origins  of  LOKANGL  can  be  traced  to  an  orbit  and  ephemeris  programs,  called 
DABOS,  developed  by  Minka,  Fein,  and  Clemenz  (1).  The  major  thrust  of  DABOS 
was  the  application  of  minimum  variance  filtering  techniques  to  raw  radar 
observations.  This  task  is  currently  performed  by  the  agencies  which  furnish 
the  orbital  element  sets  which  are  used  as  input  to  LOKANGL. 

Although  LOKANGL  avoids  filtering  raw  data  in  favor  of  a  more  computationally 
efficient  technique,  much  of  its  structure  has  been  taken  from  the  earlier 
program.  Thus  Reference  1  has  considerable  material  of  relevance  to  LOKANGL. 
Evolving  by  stages,  LOKANGL  consists  of  software  modules  of  various  origins 
which  have  been  added  as  special  needs  arose.  However,  the  basic  program  has 
now  remained  static  for  several  years.  Special  purpose  versions  continue 
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to  be  developed  (see  sections  3.0  and  4.0),  but  such  modifications  have  not 
been  made  a  part  of  the  basic  LOKANGL  system.  It  is  appropriate  at  this  stage, 
then,  to  discuss  LOKANGL,  both  to  document  the  up-to-date  version  of  the  basic 
program  and  to  form  the  basis  for  presenting  more  recent  programs  for  which 
LOKANGL  forms  the  foundation. 

The  basic  objective  of  LOKANGL  is  efficient  (in  terms  of  central  processor 
time)  calculation  of  satellite  ephemerides.  The  method  used  is  to  extrapolate 
mean  orbital  elements  by  means  of  power  series  expansions.  The  limited  range 
of  applicability  of  such  expansions  translates  into  a  definite  restriction  on 
the  extent  of  the  time  interval  over  which  valid  predictions  can  be  made. 
Nevertheless,  when  element  sets  are  available  for  epochs  sufficiently  close  to 
the  times  of  interest,  LOKANGL  provides  reliable  results  for  a  much  lower 
investment  in  computing  time  than  routines  that  require  solution  of  the  govern¬ 
ing  differential  equations. 

LOKANGL  can  be  furnished  one  or  more  sets  of  orbital  elements.  It  obtains 
the  derivatives  of  the  mean  elements,  needed  for  calculation,  in  one  of  several 
ways.  First  consider  the  case  in  which  only  one  element  set  is  furnished.  If 
the  derivatives  ,in  addition  to  elements,  are  supplied  by  the  issuing  agency, 
•■hese  values  are  used.  Otherwise,  the  time  derivatives  of  the  elements 
are  calculated  based  upon  internal  physical  models  for  the  Earth's  atmospheric 
density  and  geopotential.  Calculation  of  derivatives  based  on  these  models 
follows  procedures  presented  by  King-Hele  (2).  The  given  elements,  together 
with  their  calculated  derivatives,  are  then  used  to  perform  extrapolations 
either  backward  or  forward  in  time.  Extrapolated  mean  elements  are  then 
converted  to  equivalent  alternative  forms,  as  required  by  the  application. 
However,  when  given  two  or  more  element  sets,  and  ephemeris  data  is  required  at 
a  time  that  is  intermediate  to  the  epochs  of  two  of  these  element  sets,  LOKANGL 
makes  explicit  use  of  the  inherent  time-variation  between  these  two  successive 
sets  to  evaluate  derivatives  during  the  intervening  period.  This  feature 
overrides  the  internal  methods  of  calculating  derivatives  which  are  used  in  the 
absence  of  spanning  elements.  Any  number  of  consecutive  (i.e.,  time-ordered) 
element  sets  can  be  furnished  to  the  program  .  They  are  used  pairwise  to 
evaluate  time  derivatives  to  be  used  within  the  time  interval  spanned  between 
their  respective  epochs. 
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Satellite  ephemerides  can  be  output  in  a  variety  of  forms.  In  addition, 
various  related  quantities  can  also  be  computed.  These  include:  station  look 
angles  (i.e.,  topocentric  coordinates  of  satellite)  for  both  fixed  and  moving 
(e.g.,  an  aircraft)  stations;  coordinates  of  intersection  of  station-satellite 
1 ine-of-si ght  with  the  ionosphere;  angle  between  satellite-sun  line  and  satel¬ 
lite's  Earth  horizon  (defines  eclipse  conditions). 

Table  1  summarizes  major  features  of  the  program. 

o  Can  use  pairs  of  element  sets  to  compute  time  variation  of  ele¬ 
ments 

o  Has  internal  atmospheric  density/drag  and  geopotential  models  for 
evaluating  time  derivatives  of  elements  based  on  method  of 
Ki ng-Hel e 

o  Accepts  variety  of  types  of  orbital  elements 

o  Provides  look  angles  to  satellite  for  both  fixed  and  moving  (i.e., 
aircraft)  stations 

o  Provides  ionospheric  intersections  for  satellite-station  paths 

o  Indicates  solar  ill  umination/shadowing  of  satellite 

o  Provides  option  for  standard  output  ephemeris  file  and/or  several 
types  of  plot  files 

o  Provides  option  for  printing  variety  of  output  listings 

Table  1.  Major  Features  of  LOKANGL 


2.1.2  Approach  and  Program  Organization 


Figure  1  illustrates  the  flow  of  information  and  key  operations  in  LOKANGL. 
Table  2  exhibits  the  roles  of  key  routines.  Any  one  of  five  different  standard 
types  of  element  sets  can  be  input.  (The  type  of  elements  furnished  varies 
among  the  agencies  that  supply  them).  Element  derivatives  are  obtained 
differently  for  the  different  types  of  sets,  as  shown  in  Figure  1.  The 
internal  atmospheric  density/drag  model  and  geopotential  model  provide  data 
which  is  needed  to  perform  transformations  between  types  of  elements  and  to 
evaluate  element  derivatives.  The  interpolation  method  of  evaluating  deri¬ 
vatives  is  employed  wherever  spanning  element  sets  are  available.  Propagation 
of  elements  from  their  epoch  to  the  time  of  interest  is  performed  in  terms  of 
mean  elements.  These  mean  elements  are  then  converted  to  osculating  elements 
for  calculations  of  instantaneous  values  of  observable  parameters  such  as 
position  and  velocity.  Ancillary  calculations  include:  solar  ephemeris; 
Greenwich  sidereal  time;  ECI  to  topocentric  coordinate  transformations;  and 
solution  for  station-satellite  1 ine-of-sight  intersection  with  the  ionosphere. 

2.1.3  Input 

Input  to  LOKANGL  is  by  means  of  punched  cards.  Their  layout  and  organization 
are  shown  in  Table  3.  Note  that  five  different  types  of  element  sets  can  be 
accommodated.  Organizations  which  furnish  elements  include  the  following: 

NORAD  (mean  Keplerian  elements) 

SCF  (position/velocity  vectors) 

NASA  (Brouwer  mean  Keplerian  elements 

and  osculating  Keplerian  elements) 

COMSAT  (mean  Keplerian  elements). 

Figures  2  to  5  illustrate  the  format  of  cards  for  different  element  types. 
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Figure  1.  Simplified  Operational  and  Data  Flow  for  LOKANGL 


o  LOKANGL  (Main  Program) 


Reads  input  cards  (except  flight  cards) 

Depending  on  type  of  element  set,  calls  appropriate 
routine  to  convert  to  mean  elements  and  evaluate 
derivatives. 

If  second  element  set  is  given,  interpolates  derivatives  over 
spanning  interval 

Computes  elements  at  times  of  interest 
Calls  FLTRANS  to  read  flight  cards 

Calls  SPPROU  to  perform  ephemeris  calculations  at  each  time  of  interest 
Calls  WRSTP  to  perform  solar/ionosphere/aircraft  (moving  station)  related 
calculations 

Clocks  operation  in  incremental  steps  from  initial  to  final  time 
Terminates  job  when  final  time  is  reached. 

o  SPPROU  (Ephemeris  Routine) 

Called  for  each  type  of  calcuation 
Calculates  ephemeris  quantities 
Writes  ephemeris  on  file  T APE3 
Writes  perigee-apogee  data  on  file  TAPE7 

Converts  mean  elements  to  osculating  values  and  to  P/V  vectors 

Determines  station  viewing  parameters 

Evaluates  Greenwich  sidereal  time 

Prints  variety  of  types  of  ephemeris  output. 

o  WRSTP  (Moving  Station  Routine) 

Cal  led  only  once 

Reads  TAP E3  to  obtain  ephemeris  data 
Calls  CORFL  for  interpolated  A/C  coordinates 
Calls  SOL V I L  to  evaluate  solar  illumination  conditions 
Calls  SILL2  to  determine  station-satellite  line-of-sight  intersection 
with  ionosphere 

Prints  output  for  moving  station,  inospheric  intersection,  and 
satellite  eclipsing  calculations. 


Table  2.  Major  Routines  and  Their  Chief  Functions 


—  INPUT  FORMAT  FOR  PROGRAM  LOfCANGL  * 

-  DATA  DECK  SETUP 


CARO  COLS  OESC  SIPTION  - 

i  1  CDOE  OF  ORBITAL  DETERMINATION  FORM 

-  1=  FORM  NO.  J  SCF  2-CARJ  POS.-VEL.  SET 

2s  FORM  NO.  2  \  AOC  2-CAR}  ELEMENT  DATA  SET 

3*  3-CARD  ELEMENT  OATA  SET  - - — 

<*=OSCULATING  ELEMENTS 

S*Fu«rt  NO.  1\  AOC  5 -CARO  OATA  SET -  - 


2* 


ONE  OR  MORE  ELEMENT  SETS - - - 

8-9  X  OR  U  ON  FIRST  CARD  OF  Al.  SUBSEQUENT  ELEMENT  SETS 

a  T  INDICATES  THrUST  TIME  CARDS  ONE  CARO  SET  - 

PROVIDES  TIME  OF  THRUST  14  STANDARD  COLUMNS 


3  1-3  NO.  OF  STATIONS.  IF  NO  STATIONS  USE  0 

-----  FOR  AIRCRAFT  FLIGHT  SIMULATION  RUN,  USE  -1 - 

6  CJOE  1  OR  0  FOR  PRINT  CONTROL  OF  STATION  OATA 

0=PKINT  BT  TIME  ONLV  -  -  - 

IMPRINT  0T  STATIONS 

2-10  0  =  S  T ANOARO  - 

IONOHT  =  SUB-IONOSPH£RIC  ALT T  T UDE  (KM)  13 

- 13-15  MI  NEL  V=  HIN  ELEVATION  FOR  J  A  T  VIENING — 13 - 

16-30  FOR  AIRCRAFT  FLIGHT  SIMULATION  RUN,  ALTITUDE  OF 
AIRCRAFT  <M)  Fla.l -  -  - - 


-3* 


STATION  LOCATION  CARDS, IF  NNS.GT.O  — - - 

1-5  I.O.  OF  STATION  (NUMBER) 

—  2-8  2=GE00eriC  STSTEH  - - - 

9-23  STATION  GEOOET1C  LATITUDE.  DEGREES 

Zh-36  STATION  LONGITUOe  (POSITIV:  NEST).  DEGREES  — 

39-53  STATION  HEIGHT.  METERS 

61-22  NAME  OF  STATICN  -  - - - 


-N - 1-15  TIME  INCREMENT  IN  SECONDS.  F15.  5 - - 

12-69  TIME  INTERVAL  F CR  PRINT  OR  TAPE  3  WRITE 

-  START  COLS\  1 2- 18520-21 S23-24526-29531-34 5 36 -G1 

FINAL  CQLS\  <«3-  <.<*5<*6-<»Z;  49-50  552-555  52-605  62-62 

- -  f  I  ME  MONTHS  OAT?  TEARS  HOURS  MINS  SEC  - 

FORMAT  125  125  125  F3. 15  FJ.15  F5.3 

- FINAL  TIME  NOT  NECtSSARY  I*  NHS.EQ.-l - 


5 

5A 

3* 


AIRCRAFT  POSITION-TIME  CAROS,  IF  NHS.EQ.-l  - 

TITLE  REQUIRED  UA10I 

TWO  OR  MORE  P  CS I TI  ON-TIME  SARDS - 

I- 5  HOURS  FROM  STARTING  TIME  15 

6-10  MINUTES  F5.0  --  -  - - - - 

II- 15  OEGRSES  LATITUDE  15 
16-20  MINUTES  LATiTUOE  15 
21-25  OEGRtES  LONGITUOE  15 

25-30  MINUTES  LONGITUDE  15 

2_5  — “lOVO"  INDICATE S  ENO  OF  POS I T ION-TIME  CAROS - 


—  COOE  FOR  PRINT-OUT  0=  NO  *RINT!  1*  PRINT  FOLLOWING' 
2  IMPOSITION  ANO  VELOCITY 

4- - l  =  SUS-SATELLITc  OATA - 

6  1=  MEAN  ELEMENTS 

8  1  =  0  BSEFVAT IO  N  OATA  - - - 

•,-lU  OmSTANOARO  BINARY  ON  TAPEJ 
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11-13  1.0=SINGLE  PASS  REFRAC.  CORRECTION  1ST  AN CARD) 

-  •’ENO  OF  PROBLEM” - - - — - ~ — 


-  ABOVi  MAY  BE  REPEATED  N  TIMES 
LAST  CAROS  ”9”  IN  COL.  1 


Table  3.  Input  Format  for  Program  LOKANGL 


Figure  2.  NORAD  SPADATS  Element  Set  and  SCF  Pos i ti on/Vel oci ty  Vectors 
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Thrust  cards  are  a  special  class  of  element  card  which  signal  the  occurrence  r: 
an  orbit  adjust  at  specified  times.  The  hallmark  of  thrust  cards  ic  "he  :  m 
column  8.  Other  thrust  information  on  the  card  is: 

o  Columns  19-20,  year  ir.  format  12; 
o  Columns  21-23,  day  number  in  format  12; 
o  Columns  24-33,  time  in  seconds  in  format  F10.3. 

A  thrust  card  should  be  both  preceded  and  followed  by  standard  element  cards. 

LOKANGL  interprets  the  occurrence  of  a  thrust  as  a  discontinuity  in  the  orbit. 
The  program  divides  the  interval  between  the  epochs  of  the  element  sets  occur¬ 
ring  immediately  before  and  after  a  thrust  card  into  two  intervals: 

interval  #1  -  fr«.m  epoch  of  preceding  set  up  to  thrust  time, 

interval  #2  -  from  thrust  time  up  to  epoch  of  following  set. 

In  the  first  interval  the  program  performs  a  forward  extrapolation;  in  the 
second,  a  backward  extrapolation.  Thus,  during  both  interval  1  and  2,  it  is 
not  possible  to  obtain  derivatives  by  the  process  of  interpolating  between 
spanning  element  sets. 

2.1.4  Printed  Output 


Eight  standard  types  of  listed  output  are  available,  seven  of  which  are  select¬ 
able  by  choice  of  options  on  the  input  cards.  The  eiahth  is  a  standard  header 
which  precedes  all  listed  output.  These  output  types  are  summarized  below,  and 
examples  are  presented  in  figures  b  to  12. 

1)  Header  (common  to  all  output  print  options)  Moure  b 


Type  of  elements  used,  epoch 

Element  values  and  derivatives 

Initial  and  final  print  times;  time  increment 

Print  option  selected 

Initial  orbital  paramef  t 


2)  Sub-satellite  position  option  (0,1, 0,0  option) 


Figure  7 


Date,  time,  rev  number 

Geocentric  and  geodetic  latitude,  W.  longitude,  altitude 
Geocentric  radius,  velocity,  local  time 

3)  Mean  elements  option  (0,0, 1,0  option)  Figure  8 

Date,  time 

Semi -major  axis,  eccentricity,  i noli  nation 
Ascending  node,  argument  of  perigee,  mean  anomaly 

4)  Position/Velocity  option  (1,0, 0,0  option)  Figure  9 

Date,  time 

X,  Y,  Z,  VX,  VY,  VZ 

5)  Station  observation  option  (0,0, 0,1  option)  Figure  10 

Date,  time  (including  seconds  of  day) 

Rev  number,  elevation,  azimuth 
Range,  right  ascension,  declination 
Time  derivatives  of:  elevation,  azimuth, 
range,  declination,  right  ascension 

6)  Aircraft  flight/ionosphere  option  (NMS=-1)  Figure  11 

Flight  segment,  date,  time 

Aircraft  latitude  and  longitude 

Aircraft  to  satellite  viewing  parameters: 

elevation,  azimuth,  range,  range  rate 
Sub-ionospheric  point:  latitude,  longitude 
Sub-satellite  point:  latitude,  longitude 
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Figure  12a 


7a)  Station  look  angles/sub-satellite  point/satellite  occultation 
option  (default  option:  columns  2,  4,  6,  and  8  of  card  6 
are  all  zero;  NMSM;  and  ionospheric  height  omitted  on 
card  3) 

Station  number,  date,  time 

Station  viewing  parameters:  elevation,  azimuth,  slant  range, 
right  ascension,  declination,  altitude 

Sub-satellite  point:  latitude,  longtitude,  solar  elevation 

Illumination  (eclipse)  angle  (angle  between  satellite  -  Sun  line 
and  satellite's  Earth  horizon.) 

Rev  number 

7b)  Station  look  angles/sub-satellite  point/satellite  occultation/  Figure  12b 
sub-ionospheric  point  option  (default  option:  columns  2,  4, 

6,  and  8  of  card  6  are  all  zero;  NMS_>1;  and  non-zero  iono¬ 
spheric  height  is  entered  on  card  3) 


Station  number,  date,  time 

Station  viewing  parameters:  elevation,  azimuth,  slant  range, 
altitude 

Sub-ionospheric  point:  latitude,  longtitude 
Sub-satellite  point:  latitude,  longtitude,  solar  elevation 
Illumination  (eclipse)  angle  (angle  between  satellite  -  Sun  line 
and  satellite's  Earth  horizon.) 

Rev  number 
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Figure  7.  Sample  Sub-Satellite  Printout 
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Station  Look  Anqle/Sub-Ionosphnre/Suu-SatPll i te/Uccula tion  Printout 


If  a  user  requires  only  the  type  of  output  illustrated  by  Figure  12a  or  12b, 
columns  2,  4,  6,  and  8  of  card  6  should  all  be  zero,  and  the  deck  should 
include  at  least  one  station  card.  However,  if  the  user  elects  another  print 
option  (i.e.,  columns  2,  4,  6,  and  8  of  card  C  not  all  zero),  the  Figure  12 
printout  will  also  be  furnished  {in  addition  to  type  of  printout  specifically 
requested)  if  a  station  card  is  present  in  the  ueck. 


The  satellite  illumination  angle  shown  in  Figures  12a  and  12b  measures  the 
elevation  of  the  center  of  the  solar  disk  above  the  satellite's  Earth-horizon. 
Figure  13  illustrates  the  geometry.  The  plane  of  the  drawing  is  that  plane 
which  passes  through  the  center  of  the  Earth,  the  center  of  the  sun,  and  the 
satellite's  location.  Positive  values  of  illumination  angle  imply  that  the 
satellite  is  illuminated  by  the  Sun;  for  negative  values,  the  satellite  is 
immersed  in  the  Earth's  shadow. 


The  solar  elevation  angle  listed  on  Figures  12a  and  12b  is  simply  the  elevation 
of  the  Sun  at  the  location  of  the  specified  station. 

The  station  parameters  right  ascension  and  declination  shown  on  Figure  12a  are 
topocentric  coordinates  of  the  satellite.  They  represent  the  orientation  of 
the  station-to-satel 1 i te  vector  referenced  to  the  equatorial  plane  and  the 
direction  of  the  vernal  equinox. 


Output  Files 


The  standard  ephemeris  file  available  from  LOKANGL  is  designated  TAPE3.  It  is 
written  in  Subroutine  SPPROU.  The  format  of  this  file  is  shown  in  Table  4. 

Not  only  is  it  available  as  an  external  output  of  the  program,  but  it  is  also 
used  internally  by  Subroutine  WRS1P  to  obtain  the  parameters  necessary  to 
calculate  and  print  certain  types  of  output. 


Graphical  display  and  analysis  of  apogee/perigee  data  is  a  sensitive  tool  for 
assessing  the  caliber  of  ephemeris  calculations.  TAPE7  is  a  file  which  pro¬ 
vides  the  input  data  to  Program  PLOTIT  for  display  of  apogee/perigee  infor¬ 
mation.  The  structure  of  a  TAPE7  record  is  shown  in  Table  5. 


\ 
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Satellite 


Figure  13.  Geometry  of  tre  Illumination  Angle 
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TA^El  BINARY  FORNAT  FROOUCED  BY  SUBROUTINE  SPPROU  IS  AS  FOLLOWS 
FIRST  RECORD  IDENTIFICATION  RtCORO 

WORD  DEFINITION 


1.  <-  7 

- 2.  NJSAT - 

3.  rmsEctu 

- A.  OJUPRT(l) 

5.  TI MSEC (2 ) 

- 6,-0JUPRr<2> 

7.  OPrINT 
- a.  NDSPRI  -- 


7  i  NUMBER  OF  WOROS  REMAINING  IN  THIS  RECORD 

-SATELLITE  NUMBER— - - - - — 

TINE  OF  OAT  OF  INITIAL  PRINT  TINE  (SEC) 
-NOCIFIED  JULIAN  DATE  JF  FINAL  PRINT  TIME  — 
TINE  OF  OAT  OF  FINAL  »RInT  TINE 

-MODIFIED  JULIAN  DATE  JF- FINAt-PRIN  T -TIME - 

PRINT  INTERVAL  (SEC) 

NUMBER  OF  SPECIAL  PRINT-TINES - 


- OATm-RECORDS  — • — - — THIS  RECORO  IS  REPEATEJ -FOR  EVERT  PRINT  TIME - 


•NOR  D - DEFINITION 


n _  w  » 

riA  V  1 1  _ 

c  • 

3. 

ti**  I  JL  — ’ 

kmoout 

_ U  ;  \  A  ("N  1  l  T _ 

c 

5  • 

KY  ROUT 

- 6. 

-KHROOT 

7. 

K-IOUT 

d « 

sicour 

9. 

T'.URMO 

0/(1)  - 

li. 

0/12) 

-L|2. 

OV  <3  )  - 

13. 

OV  (A) 

-  IN. 

DV<5) 

15. 

OV  (6 ) 

— — IB. 

AL  TIO 

17. 

R«  OPRT 

-  la. 

VTOTP* 

la. 

■j  :  OC  En 

20. 

GlODET 

21. 

OLAMO 

22. 

MS  Th* 

23. 

MS  THIN 

2*. 

STSICo 

-  25. 

A*  SEMI 

2d. 

ICCEN 

27. 

XOlNCu 

2a. 

XWASC 

2  N. 

XWPERI 

II . 

XX  MEAN 

-  31. 

IKtV 

32. 

IJ 

33. 

IT  YP£ 

3*. 

NUMST  A 

35. 

cL RATE 

3d. 

AZRATE 

37. 

RA  kA  Tc 

3  o. 

DC  RATE 

3  a. 

EL  EVAT 

AC. 

azimut 

A  1. 

RANGES 

A2. 

ranrat 

AO. 

-■tasc 

A  a  • 

Cl  I  n 

The 

SEQUcnCi 

- - — NUMBER  of  horos  renainins-in  THIS-R£CORO= - 

32  ♦  11*IJ 

- NOOIFIEO  JULIAN  DATE  OF-THE-EPHEMERIS  PRINT  TIME 

CALENDAR  MONTH 

- CALENDAR  oat - 

CALENDAR  YEAR  (LAST  2  JISITS  OF  19XX) 

- —  HOUR  OF  DAT - - - —  - - - - - - - - 

MINUTE  OF  HOUR 

- —  SECONOS  OF  MINUTE - - - - - — - — - 

TIME  OF  DAY  IN  SECONDS  CORRESPONDING  TO  DAY JL 

- - X  COORDINATE  OF  POSITIJV  VECTOR  CKH) - 

T  COORDINATE  OF  POSITION  VECTOR  (KM) 

- . - 1  COO* 01  NATE  OF  POSITION  VECTOR  (KM)  - _____ 

X-DOT  COORDINATE  OF  VI.JCITY  VECTOR  (Kd/SEC) 

-  r-DOT  COORDINATE  OF  VI.OCITY  VECTOR  (KN/SEO 

Z-QOT  COORDINATE  OF  VI.OCITY  VECTOR  (KN/SEC) 

- SATELLITE  ALTITUOE  (KM) - 

SATELLITE  GEOCENTRIC  DISTRANCE  (KN) 

- VELOCITY  (KM/SEC)  - - - 

GEOCENTRIC  LATITUDE  (JIG) 


GEODETIC  LATITUOE  (DEC) 

SATELLITE  LONGITUDE  (JIG)  -  -  -  — 

HOUR  OF  GREENWICH  MEAN  SIDEREAL  TIME 

-  MINUTE  OF  GREENWICH  NUN  SIDEREAL  TINE  - 

SECONDS  OF  GREENWICH  MIAN  SICEnEAL  TINE 

—  SEMI-MAJOR  AXIS  <KN)  -  -  - - 

ECCENTRICirr 

INCLINATION  (DEG)  . 

RIGHT  ASC-NSiON  OF  ASCENDING  NODE  (DEG) 

ARGUMENT  OF  PERIGcE  (01 0) 
mean  anomaly  (ceg) 

revulut ion  NUMeER  —  - 

NUMBER  of  STATIONS  IN  MIS  DATA  RECORD  (<2l) 

*1*  FINAL  WORD  iF  I J  =  J 
STATION  number 
Elevation  ratc  (oeg/se:» 

AZIMUTH  *AT2  (DEG/CEC) 

MIGHT  ASC-nSIOn  RATc  OI6/SEC) 

OECLINATICN  RATt.  (DEG/ S  £  C) 
tLtVATION  ( OEG) 

AZIMUTH  (OEG) 

RANGE  (KM) 

RANGE  RAT-  (KM/SEC) 

RIGHT  ASCENSION  ( OEG ) 

OECLINATICN  (OEG) 

OF  HOROS  3*.  TO  a*.  IS  rEPlATI)  IJ  TIMES 


- FINmL  RECORD  - IDENT IF ICA U ON  REC0*0  -)R  THE  END  OF  PROBLEM 

—  WJRO  -  -  DEFINITION 

—  1.  KFi  1  -  - 

2.  blank  Hollerith  blanks 

An  END  OF  FILE  FJlLOWS  THIS  RlCORO 


Table  4.  Structure  of  TAPE3  Ephemeris  File 


WORD  #  QUANTITY  FORMAT 

1  Day  Number  13 

2  Hour  12 

3  Minute  12 

4  Year  12 

5  Perigee  Latitude  F6.2 

6  Perigee  W.  Longitude  F6.2 

7  Altitude  of  Perigee  (Kms)  F8.1 

8  Apogee  Latitude  F6.2 

9  Apogee  W.  Longitude  F6.2 

10  Altitude  of  Apogee  (Kms)  F8.1 

11  Rev  Number  19 

12  Local  Time  of  Perigee  (Hours)  F6.2 

13  Local  time  at  Apogee  (Hours)  F6.2 


Table  5.  Structure  of  Record  for  TAPE7,  Apogee/Perigee  File 
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Several  plotting  capabilities  have  been  developed  for  displaying  conditions  of 
eclipsing  of  the  satellite  by  the  Earth  in  support  of  the  SCATHA  Program. 

These  displays  make  use  of  the  satellite  illumination  angle,  a  quantity  that  is 
not  available  from  either  TAPE3  or  TAPE7.  It.  is  calculated  in  subroutine 
WRSTP,  which,  using  the  CDC  UPDATE  cards  shown  in  Table  6,  can  be  modified  to 
produce  TAPE2  which  contains  the  data  necessary  for  the  eclipse  plots.  Table  7 
summarizes  the  structure  of  a  TAPE2  record. 


ATTACH,  OLDPL ,  L0KPLX3421,  ID=ROBERTS ,  MR=1. 

REQUEST,  TAPE2,  *PF. 

UPDATE  (F) 

ETN  I . 

LDSET  (PRESET=ZERO) 
l GO  (Pl  =  77777  ) 

CATALOG,  TAPE2,  S PR 79 ,  ID=BREHM. 

7/8/9 
*ID  N0V6 
*D  L0KANG.3 

*TAPE2,  TAPE 7 )  (Continuation  Card) 

*1  WRSTP. 154 

367  WRITE  (2)  LMONTH,  LDAV ,  LYEAR,  JHOUR,  JMIN,  JSEC,  SLAT,  SLON,  XILLM , 
*ALTI0,  IDAY,  RADPRT,  KK  (Continuation  Card) 

7/8/9 

-  Data  Cards  - 
6/7/S/9 

-  The  following  UPDATE  modifications  can  be  used  to  suppress  listed  output  - 
*D  WRSTP -151 ,WRST* 152 
*D  WRSTP. 41 

Go  to  300 


Table  6.  Control  Cards  and  UPDATE  Cards  for  Creating  TAPE2 


r 


Word 

Number 

Variabl e 

Name 

Variable 

1 

LMONTH 

Calendar  Month 

2 

LDAY 

Calendar  Day 

3 

LYEAR 

Calendar  Year 

4 

JHOUR 

Hour  of  Day 

5 

JMIN 

Minute  of  Hour 

6 

JSEC 

Second  of  Minute 

7 

SLAT 

Sub-Satellite  Latitude 

8 

SLON 

Sub-Satellite  Longitude  (West) 

9 

XILLM 

Satellite  Illumination  Angle 

10 

ALTIO 

Satellite  Altitude  (Kms) 

11 

IDAY 

Julian  Day 

12 

RADPRT 

Satellite  Geocentric  Distance  (Kms) 

13 

KK 

Station  Number 

Table  7.  Structure  of  TAPE2  Record 
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2.1.6  Error  Checks 


Several  safeguards  have  been  built  into  the  program  to  diagnose  errors  commonly 
encountered  in  practice.  They  are  listed  in  Table  8. 


Error  Condition 

Checks 

Internal  and  External 

Indications 

Recovery  Procedures 

or  Action  Required 

Thrust  Cards 

"Error-Thrust  Cards 

out  of  sequence" 

Correct  Input  Data  Cards 

Number  of 

Stations 

"Station  printout 

requested,  NMS 
must  not  equal  zero" 

Correct  Input  Data  Cards 

Time  Interval 

"Print  interval  is 

negative  or  zero" 

Correct  Input  Data  Cards 

Print  Times 

"Start  print  time  is 

greater  than  end  print 

time" 

Correct  Input  Data  Cards 

NOSPRI 

"NOSPRI,  number  of 

special  points  is 

negative  or  exceeds  120" 

Correct  Input  Data  Cards 

ID  of  Stations 

on  TAPE3  correspond 
with  input  ID's? 

"Station  which  is  on  TAPE3 

does  not  agree  with  any  of 
the  INPUT  stations.  Program 

Identify  source  of 

erroneous  ID  on  TAPE3 

and  correct 

cannot  continue 


Table  8.  Error  Checks 
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There  are  no  graphics  within  LOKANGL  itself.  However  TAPE2  and  TAPE7  are 
intended  specifically  for  plotting  applications.  TAPE2  is  used  for  those 
plotting  applications  that  involve  evaluation  of  elipsing  conditions.  TAPE 3, 
the  general  purpose  ephemeris  file,  is  also  useful  for  various  plotting  pur¬ 
poses. 

TAPE7  is  used  in  conjunction  with  Program  PLOTIT  to  generate  apogee/perigee 
plots.  These  displays  can  be  used  as  a  sensitive  diagnostic  for  evaluating  the 
quality  of  ephemeris  calculations,  especially  for  those  cases  in  which  a  long 
time  interval  is  spanned  by  multiple  element  sets.  Non-physical  features  in 
these  plots,  _,uch  as  discontinuous  derivatives,  often  reveal  errors  in  the 
corresponding  element  sets. 

Figure  14  is  an  example  of  a  plot  generated  by  PLOTIT.  Table  9  illustrates  the 
setup  of  the  input  deck  for  PLOTIT.  Table  10  illustrates  a  typical  run-deck 
for  a  combined  L0KANGL/PL0TIT  run.  Compiled  versions  of  both  programs  have 
been  assumed  available. 
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TabVe  10.  Setup  of  Deck  for  Combined  L0KANGL/PL0T1T  Run 


87 


-  S  Nu( 


2.1.8  Computer  Requirements 

As  currently  configured,  LOKANGL  requires  67K  octal  core  memory.  Compilation 
uses  15  seconds  of  central  processor  time.  The  run  time  for  a  given  job  can  be 
estimated  by: 

Tfinal  “  ^initial 

Run  Time  =  _  x  (C)  seconds 

t 

where  t  is  the  time  step  size 

and  C  is  a  number  between  0.01  and  0.04  (CP  time  for  one  calculation). 

2.2  Sol ar-Magnetospheric  Coordi nates 

2.2.1  Introduction 

The  geocentric  solar  magnetospheric  (SM)  coordinate  system  is  defined  to  be  a 
set  of  rectangular  axes  with  the  X-axis  pointing  in  the  direction  from  the 
Earth's  center  to  the  Sun.  The  Y-axis  is  perpendicular  to  the  plane  containing 
the  X-axis  and  the  Earth's  magnetic  dipole  axis.  The  Z-axis  is  directed  in  the 
same  sense  as  the  nothern  magnetic  pole. 

In  this  system  the  Y-axis  is  always  in  the  magnetic  equatorial  plane  and  lies 
in  the  dawn-dusk  meridian,  oriented  towards  dusk.  As  the  Earth's  magnetic 
dipole  rotates  about  the  Earth's  rotational  axis,  this  system  of  coordinates 
rocks  about  the  solar  direction  with  a  24-hour  period.  An  important  feature  is 
that,  in  this  system,  the  otherwise  three  dimensional  motion  of  the  Earth's 
dipole  is  reduced  to  motion  in  a  plane  (the  X-Z  plane).  The  system  finds 
application  in  areas  of  magnetospheric  physics  such  as  the  interaction  of  the 
solar  wind  with  the  Earth's  magnetic  field. 

To  exploit  the  advantage  of  SM  coordinates,  a  capability  for  expressing  satel¬ 
lite  ephemerides  in  this  system  has  been  incorporated  into  a  version  of  Program 
LOKANGLE:  SMLOK . 
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2.2.2 


Satellite  position  and  velocity  in  Earth  Centered  Inertial  (ECI)  rectangular 
coordinates  are  available  within  subroutine  SPPROU  (See  Section  2.1).  For  the 
SM  modification,  a  software  module  which  transforms  from  ECI  to  SM  coordinates 
has  been  added  to  SPPROU.  The  SM  coordinates  are  then  written  to  a  modified 
ephemeris  file,  TAPE3,  which  can  then  be  made  available  for  purposes  external 
to  SMLOK . 

2.2.3  Input/Output 

Input  to  SMLOK  is  identical  to  that  of  a  normal  LOKANGL  run  (see  Section 
2.1.3).  The  output  listing  is  shown  in  Figure  15.  TAPE3,  the  SMLOK  ephemeris 
file,  has  the  structure  of  the  normal  LOKANGL  TAPE3,  except  that  ECI  rectangular 
coordinates  are  replaced  by  solar  magnetospheric ,  preserving  the  X,  Y,  Z  ordering. 

A  modified  version  of  the  SCATHA  5-in-l  plotting  routine  has  been  used  to  plot 
SM  coordinates.  An  example  is  shown  in  Figure  16  for  the  SCATHA  vehicle. 
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Figure  15.  Sample  Solar/Magnetospheric  Ephemeris  printout 


SCRTHfl  2/  4/79  OHR  36MN 


Figure  16.  Display  of  Solar/Magnetospheric  Coordinates 
for  the  SCATHA  Satellite 


2.2.4  Mathematical  Method 


The  subroutine  SOLMAG  converts  a  set  of  ECI  coordinates  and  associated  uni¬ 
versal  time  to  geocentric  solar  magnetospheric  coordinates.  The  transformation 
is  based  upon  four  quantities:  Greenwich  sidereal  time,  solar  right  ascension 
and  declination,  and  the  orientation  of  the  Earth's  magnetic  north  pole. 

Greenwich  sidereal  time  is  calculated  using  an  algorithm  for  ephemeris  sidereal 
time  from  Reference  (3).  The  time  is  then  converted  from  the  ^phemeris 
meridian  to  the  Greenwich  meridian  by  a  linear  transformation.  The  solar  right 
ascension  ana  declination  are  calculated  by  means  of  routines  g\  en  in  References 
(3)  and  (4). 


These  parameters  are  combined  with  the  position  of  the  Earth's  magnetic  north 
pole  in  the  following  set  of  vector  equations  defining  solar  magnetospheric 
coordinates  (AX,  AY,  AZ). 


Let  THET  and  PHI  define  the  latitude  and  longitude  of  geomagnetic  north  pole 
and  GST  be  Greenwich  sidereal  time.  Then  the  Earth  centered  inertial  co¬ 
ordinates  (ECI)  of  the  unit  vector  having  the  orientation  of  the  geomagnetic 
north  pole  are 

(DX\  /COS(THET)  COS  (GST  +  PHI) 

DY  )  =  (  COS(THET)  SIN  (GST  +  PHI ) 

DZ/  \  SIN(THET) 


Also,  let  the  position  of  the  Sun's  declination  and  right  ascension  be  DECS  and 
RAS.  Then  the  ECI  coordinates  of  the  unit  vector  oriented  in  the  direction  of 
the  Sun  is 


S= 


(COS(DECS)  COS(RAS) 
COS(DECS)  SIN(RAS) 
SIN(DECS) 
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Now  define  a  vector  T  = 


as 


ft  =  D1  x  St, 


where  yt  denotes  the  transpose  of  T. 


Also  define  U  as 


Ut  *  S11  xT1. 

Then  if  (AS,  AT,  AU)  denotes  ECI  coordinates,  the 
coordinates  are: 


solar/magnetospheric 
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3.0 


Ionospheric  Research  Support 

In  the  early  history  of  ionospheric  research,  observations  were  generally 
limited  to  those  made  from  fixed  ground  stations.  This  restricted  researchers' 
ability  to  track  specific  moving  features  with  passing  time  and  to  view  struc¬ 
tures  from  different  aspects.  Today,  satellites  and  instrumented  aircraft 
constitute  moving  platforms  which  are  used  to  enable  researchers  to  perform  a 
broader  range  of  experiments.  Time  correlated  observations  are  made  simul¬ 
taneously  from  multiple  observing  platforms.  In  planning  for  these  experi¬ 
ments,  and  analyzing  the  data,  it  becomes  necessary  to  evaluate  a  variety  of 
viewing  geometries,  especially  as  they  change  with  time.  Two  programs,  SKY  and 
IONTRK,  have  been  developed  to  satisfy  this  type  of  requirement. 

3.1  Program  SKY 

3.1.1  Introduction 

A  very  fruitful  technique  of  experimental  ionospheric  research  is  the  recording 
of  optical  radiation  associated  with  ionization  present  in  the  upper  atmosphere. 
Optical  systems  are  employed  that  provide  360°  of  azimuthal  coverage  together 
with  broad  elevation  coverage  that  extends  downward  from  the  zenith. 

Program  SKY  was  written  to  display,  in  plot  format,  specific  contours  within 
the  field  of  view  of  an  observer  using  an  "all -sky"  camera  or  TV  system. 

The  camera,  mounted  upon  an  aircraft,  takes  a  series  of  photographs  of  the  sky 
covering  360°  of  azimuth.  The  computed  plots  are  superimposed  upon  these 
photographs  to  indicate,  as  benchmarks,  the  relative  locations  and  motions  of 
optical  features  appearing  within  the  field  of  view  of  the  optical  system. 

Specifically,  geographic  or  geomagnetic  latitude  and  longitude  contours  at 
particular  altitudes  are  plotted  from  a  particular  aircraft  position.  There 
is  also  the  capability  to  plot  specified  isolated  points  in  space. 
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Much  of  the  program  is  devoted  to  making  available  the  various  input  options. 
Among  these  are: 


1)  Isolated  points  or  entire  contours  can  be  plotted. 

2)  Contours  of  constant  geographic  or  geomagnetic  latitude  or  longitude 
can  be  specified. 

3)  Aircraft  position  can  be  determined  from  standard  flight  track  cards. 
Plots  are  produced  every  ten  minutes.  Typically,  two  latitude  and 
two  longitude  contours  are  projected  onto  the  field  of  view. 

4)  Aircraft  position  can  also  be  determined  from  a  physical  tape  obtained 
from  the  aircraft's  inertial  navigation  system.  The  program  operation 
in  this  case  is  the  same  as  3)  above. 

5)  Manual  input  is  available  in  which  aircraft  position  and  contours  or 
points  are  input  directly  to  the  program. 

3.2.2  Program  Function  and  Organization 

The  geometry  of  the  problem  to  be  solved  is  shown  in  Figure  1.  Table  1  lists 
the  major  functions  which  must  be  performed  by  Sk i .  The  contours  to  be  plotted 
could  be,  for  example,  a  particular  geographic  or  geomagnetic  meridian  at  a 
given  altitude.  If  the  specification  is  in  terms  of  geographic  coordinates, 
one  immediately  has  the  points  to  be  displayed  represented  by  their  latitude, 
longitude,  and  altitude.  If  the  specification  is  for  magnetic  contours  the 
points  of  interest  must  be  converted  to  the  geographic  system.  Subroutine 
C0RR6M2  perforins  this  function. 

Next,  subroutine  AZEL  is  used  to  transform  geographic  coordinates  to  equivalent 
azimuth  and  elevation  referenced  to  the  aircraft's  location.  An  instrument 
function  is  then  applied  to  transform  azimuth  and  elevation  to  the  coordinates 
of  the  display  of  the  particular  device  being  used  (all-sky  camera  or  TV 
camera).  Account  is  taken  of  the  heading  of  the  aircraft  to  insure  that 
azimuth  is  correctly  referenced  to  the  direction  of  north.  Table  1  summarizes 
the  key  functions  performed  by  SKY. 
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Figure  1.  Viewing  Geometry 


o  Input  of  A/C  position/time  history  and  specification  of  coordinates 
to  be  displayed. 

o  Conversion  of  geomagnetic  coordinates  (if  necessary)  to  geographic 

o  Conversion  of  geographic  coordinates  of  points  to  be  plotted  to 
azimuth/elevation  referenced  to  aircraft  location 

o  Plotting  of  overlays 

Table  1.  Principal  Functions  Performed  by  SKY 
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The  preceding  operations  are  performed  for  a  given  instant  of  time,  for  which 
the  location  (latitude,  longitude,  altitude)  of  the  aircraft  is  known.  The 
output  is  a  plotted  overlay.  Time  is  then  incremented,  yielding  a  new  aircracft 
location.  The  process  is  repeated,  producing  another  overlay.  When  the 
aircraft  track  is  defined  by  standard  flight  cards,  subroutines  FLTRANS  and 
CORFL,  modified  versions  of  similar  routines  used  in  Program  LOKANGL,  process 
the  flight  card  information  to  yield  latitude  and  longitude  at  the  desired 
equi-spaced  time  intervals.  These  routines  interpolate  aircraft  geographic 
positions  for  times/locations  intermediate  to  those  specified  on  successive 
flight  cards. 

Figure  2  illustrates  the  flow  of  information  in  Program  SKY  as  the  foregoing 
operations  are  implemented.  To  minimize  core  memory  requirements,  the  program 
is  organized  into  five  overlays.  Communication  between  overlays  is  by  means  of 
labelled  common  blocks  and  files  TAPE3  and  TAPE4.  The  files  contain  points  to 
be  plotted.  The  labelled  common  contains  input  and  option  selection  data.  The 
division  of  functions  among  the  overlays  is  shown  in  Figure  2,  where  dashed 
outlines  are  used  indicate  individual  overlays. 


Figure  2.  Information  Flow  Diagram  for  SKY 
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3.1.3  Mathematical  Procedures 

SKY  requires  the  transformation  between  the  coordinates  of  a  point  in  space  and 
its  image  in  the  all -sky  camera  photograph.  Subroutine  AZEL  calculates  the 
azimuth  and  elevation  of  an  observation  point  from  a  given  aircraft  position. 
The  subroutine  uses  the  following  method  to  calculate  elevation,  EL,  and 
azimuth,  AZ. 


Let 


and 


X  =  (Xj,  X2,  X3)  be  aircraft  position 


Y  =  (Yj,  Y2,  Y3)  be  observation  point 


relative  to  geocentric  orthogonal  rectangular  coordinate  system. 


The  elevation  is  calculated  as: 


EL  =■ y  -  ARC  COS 


(-14} 

(  1X1  x  1Y1) 


The  azimuth  is  calculated  as  follows: 


Let 


N  =  (0,0,1)  be  unit  north  polar  vector. 


and 


N  x  X  (-X2,  Xi,  0) 

E  =  - -  =  - - - 

IXI  IXI 


be  eastward  directed  unit  vector 


1 
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Then  the  local  north  is: 


LN  =  E  x  X 
IX! 

The  azimuth  can  be  calculated  as: 


(Y  -  X)  .  E 
AZ  =  ARC  TAN  _ 

(  (Y  -  X)  .  LN 

3.1.4  Input 

Input  of  choice  of  options  is  by  means  of  cards.  Aircraft  navigation  information 

can  be  furnished  either  by  means  of  flight  cards  or  by  a  physical  tape;  the 
instrument  navigation  system  (INS)  tape.  The  set-up  of  input  cards  is  shown  in 
Table  2. 
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Card 

No. 

Vari able 

Name 

Format 

Variable  Description 

1 

IOPT 

511 

Option  Parameter  Array 

2 

RADIUS 

F10.2 

Plot  Radius 

XMINEL 

F10.2 

Horizon  Elev  Angle 

X10NHEI 

F10.2 

Observation  Height 

THESE 

CARDS  ARE  FOR  MANUAL 

MODE  OF  OPERATION. 

3 

LMON 

12 ,  IX 

Month 

LDAY 

12, IX 

Day 

LYR 

14, IX 

Year 

LHR 

12, IX 

Hour 

LMIN 

12, IX 

Mi nute 

LSEC 

12 

Seconds 

4 

ACLAT 

F10.2 

A/C  Latitude 

ACLON 

F10.2 

A/C  Longitude 

ALT 

F10.2 

A/C  Altitude 

ARHEAD 

F10.2 

A/C  Heading 

5 

I 

I1.9X 

No.  of  Contours 

PLLAT 

6F10.2 

Latitude  of  Contours 

6 

J 

11, 9X 

No.  of  Contours 

PILLON 

6F10.2 

Longitude  of  Contours 

Repeat  cards  3  to  b  for  each  successive  overlay  ending  program  with  LM0N=99  on 
card  3. 


For  A/C  position  from  cards  ( IOPT ( 1 )=2 )  insert  flight  track  cards  after  card  2 
ending  with  =1000. 


10PT  is  an  option  parameter  array  controlling  program  execution. 


I0PT(1) 
IOPT (2 ) 
IOPT (3 ) 
IOPT (4 ) 
IOPT  (5) 


1  A/C  Position  Taken  from  TAPE2 

2  A/C  Position  Taken  from  Flight  Cards 

1  Geographic  Coordinates  Plotted 

2  Geomagnetic  Coordinates  Plotted 

1  TV  Camera  Lens  (linear) 

2  All  Sky  Camera  Lens  (nonlinear) 

1  Latitude  and/or  Longitude  Contours  Plotted 

2  Isolated  Points  are  to  be  Plotted 
2  Manual  Mode  of  Operation 


Table  2.  Input  Cards  for  SKY 


103 


3.1.5  Output 


The  chief  output  of  Program  SKY  is  a  series  of  plots,  drawn  to  the  same  scale 
as  the  photographic  images  with  which  they  are  to  be  used  as  overlays.  A 
typical  example  is  shown  in  Figure  3.  Within  the  circular  frame  are  the 
contours  which  are  intended  to  be  combined,  as  an  overlay,  with  a  photograph 
from  the  "all-sky"  optical  system.  The  zenith  maps  into  the  point  at  the 
center  of  the  circle;  and  radial  distance  from  this  point  corresponds  to  zenith 
angle.  Registration  for  azimuth  angle  is  provided  by  the  labelled  axes. 

A  printout  is  provided  that  tabulates  the  track  of  the  aircraft.  As  an  error 
check,  the  desired  input  options  should  be  checked  against  those  echoed  at  the 
top  of  the  printout. 

3.1.6  Program  Restrictions  and  Timing 

Grids  larger  than  5.0  inches  will  cause  information  to  exceed  the  plotting 
dimension  in  the  Y-  direction.  200  inches  of  plot  is  the  limit  in  the  X- 
direction.  If  more  is  needed,  divide  the  input  into  two  smaller  runs.  Check 
input  echo  at  top  of  output  to  see  if  this  matches  the  user's  intentions. 

SKY  requires  110000  octal  core  memory.  A  typical  time  allocation  is  100 
seconds.  Control  card  set-up  for  a  typical  run  is  shown  in  Table  3. 
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Figure  3.  Typical  Overlay  Plot  from  Program  SKY 
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LOBOB ,  TlOO ,  CMllOOOO .  (Job  card) 

REQUEST,  PLOT,  *Q. 

ATTACH,  PEN,  ONLINE  PEN. 

LIBRARY,  PEN. 

ATTACH,  TAPE1 ,  MJDATA,  ID  =  ROBERTS,  MR=1. 
LDSET,  PRESET  =  ZERO. 

LGO. 

DISPOSE,  PLOT,  *PL . 

7/8/9  (END  OF  RECORD) 

[INPUT  CARDS] 

7/8/9  (END  OF  RECORD) 

6/7 /8/9  (END  OF  INFORMATION) 


Table  3.  Typical  Control  Card  Sequence  for  SKY 
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IONTRK 


3.2.1  Background 

The  complexity  assumed  by  some  ionospheric  experiments  is  evident  in  the  case 
in  which  two  instrumented  aircraft  make  simultaneous  observations.  The  first 
forms  one  leg  of  a  satel 1 ite-to-ai rcraft  radio  path.  Interest  centers  in 
interactions  of  the  radio  wave  with  the  ionosphere  through  which  it  passes.  A 
standard  function  of  Program  LOKANGL  is  identification  of  the  location  of  this 
interaction  region  in  terms  of  the  "sub-ionospheric  coordinates".  The  second 
aircraft  is  to  make  observations  of  this  interaction  region.  For  purposes  of 
experiment  planning,  and  for  subsequent  data  analysis,  the  experimenter  needs 
to  know  the  azimuth  and  elevation  angles,  at  the  second  aircraft,  of  the  iono¬ 
spheric  interaction  region  for  the  first  aircraft.  The  geometry  is  depicted  in 
Figure  4.  The  IONTRK  version  of  LOKANGL  was  developed  to  satisfy  this  require¬ 
ment.  The  reader  may  wish  to  refer  to  Section  2  for  background  information  on 
LOKANGL. 


3.2.2  Approach 

Necessary  input  to  the  calculation  is  the  evaluation  of  the  subionospheric 
point  associated  with  the  first  aircraft.  To  obtain  this  data,  the  program 
performs  what  is  essentially  a  normal  LOKANGL  sequence  of  operations;  however 
the  output  of  subroutine  WRSTP  is  written  to  file  TAPE1. 

Control  then  returns  to  the  main  routine,  LOKANGL,  as  usual.  It  is  at  this 
point  that  the  major  modifications  occur.  Normally  LOKANGL  would  proceed  to 
terminate.  This  chain  of  commands  is  interrupted  and  the  new  calculations  are 
inserted,  as  shown  in  Figure  5.  FLTRANS  is  called  to  read  and  process  data 
from  a  second  set  of  flight  cards  which  correspond  to  the  track  of  the  second 
aircraft.  This  information  is  written  to  TAPE2.  TAPE1  is  read,  to  make 
available  within  LOKANGL  the  data  relating  to  aircraft  one.  CORFL  is  then 
called  to  read  TAPE  2  and  provide  interpolated  latitudes,  longitudes,  and  times 
for  aircraft  two.  At  this  point,  the  coordinates  of  both  the  viewing  system 
and  the  region  being  viewed  are  defined  for  the  specified  instant  of  time. 
Calculation  of  the  azimuth  and  elevation  can  now  proceed.  The  results  are  then 
printed  out  and  the  program  terminates. 
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Normal  LOKANGL 
Sequence  to 
this  Point 


Figure  5.  Flow  Chart  for  Program  IONTRK 
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3.2.3  Mathematical  Method 


The  problem  of  evaluating  the  azimuth  and  elevation  is  that  addressed  by 
subroutine  AZEL  described  in  Section  3.1.3,  to  which  the  reader  is  referred. 

3.2.4  Input/Output 

The  organization  of  the  data  deck  for  IONTRK,  which  is  shown  in  Table  4,  is 
basically  similar  to  that  for  LOKANGL.  The  major  differences  are  the  following 
a  modification  to  card  3  to  enter  both  the  variable  IONTRK,  which  selects  the 
IONTRK  mode;  and  the  altitude  of  the  second  aircraft.  In  addition  the  card  7 
series,  which  represents  flight  cards  for  the  second  aircraft,  has  been  added. 

A  sample  output  is  shown  in  Figure  6.  The  flight  segment  referred  to  in  the 
first  two  columns  is  defined  as  that  portion  of  the  aircraft  flight  path 
intermediate  between  two  successive  flight  cards.  Segment  number  one,  for 
example,  falls  between  the  points  given  on  the  first  and  second  flight  cards. 
Looking  downward,  azimuth  is  measured  counterclockwise  from  north. 
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Figure  6.  Sample  Printout  from  Program  IONTRK 
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SCATHA  Support 
4.1  Introduction 


The  electrical  charging,  both  absolute  and  differential,  of  spacecraft  in 
orbit  has  been  found  to  have  a  potentially  adverse  effect  on  performance  of 
satellite-borne  systems.  Among  the  objectives  of  the  SCATHA  experiments  are 
the  investigation  of  the  charging  process  and  methods  for  its  control. 

Important  among  natural  phenomena  affecting  the  charging  process  is  the 
flow  of  photoelectron  currents  induced  by  incident  solar  ultraviolet  radiation. 

An  event  having  a  dramatic  effect  on  charging  is  suppression  of  illumination 
during  eclipsing  of  a  vehicle  as  it  passes  within  the  Earth's  shadow.  Another 
natural  environmental  factor  affecting  charging  is  the  interaction  between  the 
vehicle  and  the  dense,  hot  plasmas  within  which  it  may  be  imbedded  during 
magnetospheric  substorms.  Thus,  particle  environments  and  solar  illumination 
conditions  are  of  importance  in  interpreting  SCATHA  experimental  measurements. 

The  paths  of  energetic  charged  particles  in  the  magnetosphere  and  upper 
ionosphere  tend  to  be  guided  by  the  lines  of  the  Earth's  magnetic  field,  along 
which  the  particles  spiral.  (1,2)  The  field  lines,  then,  provide  conducting 
channels  for  the  currents  which  these  flowing  particles  create.  In  order  to 
compare  SCATHA  observations  with  data  from  other  sensors,  it  is  of  interest  to 
identify  time  intervals  when  the  SCATHA  vehicle  is  magnetically  coupled  to 
ground  stations  and/or  to  other  space  vehicles  of  interest. 

To  fulfill  this  need,  two  software  systems  have  been  developed.  The  first,  CONJ/ 
FPLOT,  traces  the  movement  in  the  northern  hemisphere  of  the  magnetic  footprint 
of  SCATHA.  The  precise  footprint  of  interest  here  is  the  one  at  ionospheric 
height,  usually  100  km,  where  the  highly  conductive  E-layer  connects  with 
the  conducting  path  along  the  magnetic  field  lines. 

The  second  system,  FTPRNT,  searches  for  periods  when  SCATHA  and  a  second 
vehicle  intersect  a  common  field  line.  Under  this  condition,  the  particle 
environments  of  the  two  vehicles  should  be  closely  related.  The  third  system. 
Program  FOOT,  evaluates  solar  illumination  conditions  at  footprint  locations. 


113 


The  change  in  illumination  conditions  as  SCATHA  passes  into  and  out  of  the 
Earth's  shadow  should  correlate  well  with  variation  in  vehicle  potential  at 
times  when  the  vehicle  is  immersed  in  a  dense  plasma  of  high  energy  particles 
associated  with  magnetospheric  substorms.  To  assess  this  relationship, 
calculations  at  several  levels  of  precision  were  performed.  Program  PENUMB 
provides  times  of  umbra!  and  penumbra!  passage.  LOKANGL ,  coupled  with  a 
special  plotting  package,  ECPLT,  provides  summary  plots  of  eclipse  conditions 
during  the  two  annual  seasons:  fall  and  spring. 

The  Air  Weather  Service  has  provided  a  variety  of  astro-geophysical  back¬ 
ground  data  to  support  SCATHA  analyses.  This  information  includes  particle 
density  distributions  with  respect  to  energy.  These  data  come  packed  on  tape 
prepared  by  a  UNIVAC  1110  system.  A  program  was  developed  to  unpack  this  data 
and  display  it  on  microfiche  plots. 

4.2  Program  C0NJ/FPL0T 

4.2.1  Objective 

The  C0NJ/FPL0T  syscem  is  intended  to  accept,  as  input,  an  ephemeris  tape, 
in  LOKANGL  TAPE3  format,  and  to  provide,  as  output,  a  plot  of  the  time  history 
of  movement  of  the  100km  footprint  of  the  SCATHA  satellite  in  the  northern 
hemisphere.  The  plot,  an  overlay  in  polar  form,  uses  as  coordinates  the 
corrected  geomagnetic  latitude  and  longitude.  The  scale  should  match  that 
of  an  existing  Northern  polar  projection  map  of  the  world,  depicting  continental 
outlines  in  corrected  geomagnetic  coordinates.  The  footprint  plot  is  to  be 
used  as  an  overlay  for  this  map,  and  includes  appropriate  registration  marks 
at  the  center  and  along  the  prime  magnetic  meridian. 

4.2.2  Approach 

Depending  on  the  overall  time  period  to  be  covered  and  the  time  interval 
between  successive  ephemeris  calculations,  Program  CONJ  accepts  either  one  or 
two  ephemeris  files  (TAPE  2,  TAPE  3).  Figure  1  shows  the  flow  of  operations. 
Subroutine  MGFLD2  is  called  to  perform  the  field  line  tracing.  (Actually, 
this  function  is  performed  by  LINTRA,  a  subroutine  in  the  MGFLD2  package.) 
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The  line  is  traced  both  up  and  down  the  field  line  from  the  vehicle's  location. 
It  is  not  obvious,  a  priori,  which  will  yield  the  footprint  in  the  northern 
hemisphere.  When  the  geographic  coordinates  of  the  footprint  have  been  deter¬ 
mined,  subroutine  C0RRGM2  is  called  to  perform  transformation  to  corrected 
geomagnetic  coordinates.  The  results  are  printed  in  tabular  form  and  written 
to  file  TAPE  11. 

A  continuous  plot  of  footprints  has  been  found  to  be  too  complex  to  be  useful, 
as  data  for  successive  days  overlap  extensively.  This  problem  has  been  avoided 
by  presenting  footprints  plotted  only  for  a  selected  day  of  the  week  (e.g.  for 
all  Mondays,  but  only  Mondays,  within  the  overall  time  period  of  interest.) 

The  second  phase  of  the  operation  involves  excercising  program  FPLOT  to 

read  file  TAPE11  and  plot  the  resulting  points.  CONJ/FPLOT  are  run  sequentially 

in  one  pass,  linked  together  by  job  control  card s. 

4.2.3  Input/Output 

The  set-up  of  input  cards  for  CONJ  is  shown  in  Table  1.  Typical  printed 
output  is  shown  in  Figure  2.  Input  and  output  files  for  CONJ  are  listed  in 
Table  2. 

Input  cards  for  FPLOT  are  described  in  Table  3.  FPLOT  produces  a  summary 
listing,  shown  in  Table  4,  which  echoes  the  specified  input  conditions  and 
summarizes  key  plotting  data  for  each  successive  day  plotted. 

Figure  3  is  an  example  of  the  plot  which  is  the  final  product.  Notice  that 
each  hour  is  marked  by  a  dot  and  the  latest  time  plotted  for  given  day  is 
denoted  by  a  dot  within  a  circle.  The  earliest  time  is  referenced  to  the 
first  digit  of  the  day  of  the  month  of  the  date. 
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Card  No. 

Variable 

Columns 

Format 

Description 

1 

NMON 

1-5 

15 

Beginning  Month 

NDAY 

5-10 

15 

Beginning  Day 

NYR 

11-15 

15 

Beginning  Year 

MYDAY 

16-20 

15 

Corresponding  Weekday  No 
Where  1  =  MON.  through 
7  =  SUN.  (0<MYDAY<8) 

HALT 

21-30 

F1O.0 

Geodetic  Altitude  of 
Conjugate  Intersection 

2 

COEF 

1-20 

2A10 

Coefficient  File  Label 

TM 

21-30 

F10. 2 

Update  Year  for  Coef. 
File  (19xx.xx) 

Table  1.  Format  of  Input  Cards 

for  Program  CONJ 
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Figure  2.  Printout  of  Location  of  Footprint  of  Satellite 

in  Geographic  and  Corrected  Geomagnetic  Coordinates 

Copy  available  to  DTiC  does  not 

permit  fully  legible  reproduction 


o  TAPE 
0  TAPE 
o  TAPE 

o  TAPE 

o  TAPE 


Card  No. 
1 

2 


1  IGRF 75  magnetic  field  coefficients  used  by  MGFLD  2. 

2  Vehicle  ephemeris  file  in  LOKANGL  TAPE3  format. 

3  Additional  ephemeris  file,  if  required;  sequential  with 
TAPE2. 

10  Corrected  geomagnetic  field  coefficients  used  by  C0RRGM2; 
Gustafsson  or  Hakura 

11  Output  file,  CGLAT  and  CGLON  of  footprints. 

Table  2.  Files  used  in  Program  CONJ 


Variable  Format 

ADAY  (1-16)  16F5.0 

MON  (1-8)  8A10 

Table  3.  Format  of  Input  Cards  for 


Description 

Day  of  Month  Values 
for  Weekday  (MYDAY ) 

Months  Corresponding 
to  ADAY  (Left  Justified) 

Program  FPLOT 


i 
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flue  li) 


4.2.4  Control  Cards  and  System  Requirements 

The  following  is  the  sequence  of  control  cards  used  to  run  CONJ  and  FPLOT 
sequentially  in  a  single  pass. 

LPOL ,  T300,  CM1 30000. 

REQUEST,  TAPE11,  *PF. 

FTN,  SL,  R=3 . 

ATTACH,  TAPE1 ,  IGRF75X342L,  ID=L0GBUC,  M?  =  l. 

ATTACH,  TAPE2,  SCAJUL10AUG13,  ID=L0LEW,  CY=1,  W  =  l. 

ATTACH,  TAPE3,  SCAAUG15SEP19 ,  ID=L0LEW,  CY=1 ,  MR=1. 

ATTACH,  TAPE10,  MJDATA,  1D=L0GI CON ,  MR  =  1 . 

ATTACH,  XX,  IORLIB,  ID  =  LOGICON,  MR  =  1 . 

LIBRARY,  XX. 

LDSET ,  PRESET=ZERO . 

LGO. 

REWIND,  TAPE  1 1 . 

RETURN, LGO. 

FTN,  SL,  R  =  3. 

REQUEST.  PLOT,  *Q. 

ROUTE,  PLOT,  DC  =  PL,  DEF. 

ATTACH,  PEN,  ONLINEPEN. 

LIBRARY,  PEN. 

LDSET,  PRESET=ZERO . 

LGO. 

A  typical  CONJ/FPLOT  run  consumes  300  seconds  of  central  processor  time 
and  requires  130K  octal  core  memory. 


4.3  Program  FTPRNT 


4.3.1  Objective 

FTPRNT  is  a  program  which  has  been  developed  to  identify  time  periods  during 
which  the  SCATHA  vehicle,  in  its  near-synchronous  altitude  orbit,  and  DMSP  F2, 
in  its  much  lower  altitude  orbit,  intersect,  within  a  suitable  tolerance,  a 
common  line  of  the  Earth's  magnetic  field. 

4.3.2  Approach 

A  significant  feature  of  the  problem  is  the  disparity  in  altitudes,  and  hence 
orbital  periods,  between  the  two  vehicles.  This  causes  an  asymmetry  between 
the  handling  of  the  vehicles  in  the  program.  Conjugate  footprints  are  traced 
from  the  location  of  the  higher  altitude  vehicle  down  to  the  altitude 
of  the  lower  altitude  vehicle  at  that  time.  Thus,  the  two  footprints  and  the 
lower  altitude  vehicle  lie  on  a  common  spherical  shell.  The  quality  of 
"closeness"  of  the  two  field  lines  intersecting  the  satellites  is  quantified 
in  terms  of  the  spatial  separation  on  this  shel 1  of  the  SCATHA  footprints  arc 
the  position  of  the  lower  satellite.  Intersection  is  defined  to  occur  wheneve 
the  vehicle  -  footprint  separation  in  the  shell  subtends  an  angle  at  the 
center  of  the  Earth  which  is  less  than  five  degrees  in  magnitude. 

4.3.3  Functional  Description 

The  program  is  divided  into  two  major  parts.  A  driver  program  FTPRNT  reads 
two  time-coincident  ephemerides  files  created  by  LOKANGL  in  the  TAPES  format. 
It  then  calls  MGFLD2  to  calculate  the  coordinates  of  the  conjugate  footprints 
of  SCATHA  at  the  altitude  of  the  lower  satellite. 


Next,  a  subroutine,  SEARCH,  is  called  at  successive  instants  of  time  to 
identify  when  the  lower  satellite  passes  within  b°  of  either  of  the  SCATHA 
conjugate  footprints.  In  order  to  interpolate  the  ephemeris  and  simplify  the 
analysis,  it  is  assumed  that  the  conjugate  footprints  of  SCATHA  move  slowly 
with  respect  to  the  motion  of  the  DMSP  F2  satellite  and  that  the  altitude  of 
DMSP  F2  does  not  change  appreciably  over  a  span  of  six  time  intervals  within 
the  ephemeris.  These  assumptions  can  be  satisfied  by  choosing  a  time  increment 
for  the  DMSP  F2  ephemeris  that  is  sufficiently  small  (say,  3-10  minutes). 

4.3.4  Mathematical  Method 

Subroutine  SEARCH  performs  the  principal  analysis  within  FTPRNT.  It  identi¬ 
fies  times  and  positions  when  two  satellites  intersect  a  common  magnetic 
field  line.  The  routine  is  used  within  the  driver  routine  which  reads  two 
time-coincident  ephemerides  tapes.  The  higher  altitude  satellite  is  placed  on 
file  TAPE2;  the  lower,  on  TAPE3.  The  routine  MGFLD2  is  called  tu  determine 
the  position  of  the  two  conjugate  points  of  the  higher  altitude  satellite 
traced  along  the  field  line  to  the  altitude  of  the  lower  altitude  satellite. 
This  information  is  passed  to  the  subroutine  SEARCH  which  then  determines  if 
the  lower  altitude  vehicle  passes  within  five  degrees,  measured  from  Earth's 
center,  of  either  of  the  conjugate  point  positions.  Such  events,  when  identi¬ 
fied,  are  printed  by  subroutine  SEARCH. 

The  routine  tabulates  distances  between  conjugate  points  and  the  lower  alti¬ 
tude  satellite  searching  for  minimum  arc  distance.  When  a  minimum  is  found  a 
check  is  made  to  determine  if  this  distance  is  below  a  preassigned  tolerance. 

If  this  condition  is  satisfied,  a  calculation  is  made  to  interpolate  the 
actual  minimum  between  the  tabulated  values. 

Here  two  assumptions  are  made  concerning  the  relationship  between  the  two 
satellites.  First,  it  is  assumed  that  the  conjugate  point  moves  slowly  with 
respect  to  the  lower  altitude  satellite.  In  fact,  it  is  assumed  that  the 
conjugate  point  remains  stationary  over  three  consecutive  times  of  calculation 
centered  at  the  tabulated  minimum.  Second,  it  is  assumed  that  over  a  space 
of  three  consecutive  times  the  positions  of  the  lower  altitude  satellite  are 
located  on  an  earth  centered  sphere  of  radius  the  same  as  the  tabulated 
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mi ni mum. 


These  assumptions  simplify  the  analysis  considerably  and  are  reasonable 
for  many  vehicles  provided  the  time  interval  between  tabulated  positions  is 
chosen  sufficiently  small. 

The  problem,  which  is  illustrated  in  Figure  4,  can  now  be  solved  by  means  of 
spherical  trigonometry.  The  appropriate  equations,  derived  for  the  pair  of 
right  spherical  triangles,  are: 

COS  (A)  =  COS(ARCF)  COS(ARCMIN) 

COS  (B)  =  COS(E-ARCF)  COS(ARCMIN), 

which  can  be  combined  as. 


tan (ARCF  )  = 


COS  B 
COS  A 


COSE 

sin  E 


Also,  the  minimum  arc,  ARCMIN,  can  be  expressed  as, 


COSA 

COS  (ARCMIN)  =  - 

COS  ARCF 

Now  an  interpolation  can  be  made  to  determine  the  actual  position  and  time  of 
the  minimum,  as  well  as  the  position  and  time  when  the  lower  altitude  satellite 
passes  within  five  degrees  of  the  conjugate  point. 

4.3.5  Input/Output 

Input  for  FTPRNT  consists  of  punched  cards  and  two  data  files.  The  cards 
required  and  their  format  are  shown  in  Table  5.  The  two  files,  typically 
permanent  files,  are  the  ephemerides  for  the  high  and  low  altitude  satellites, 
e.g.,  SCATHA  and  DMSP  F2.  These  are  in  the  standard  LOKANGL  TAPE3  format. 

A  sample  output  listing  is  shown  in  Figure  5. 
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Card 

No. 

Variable 

Name 

Card 

Col. 

Format 

Variable  Description 

1 

COEF (2) 

1-20 

2A10 

Coefficient  Label 

1 

TM 

21-30 

F10.2 

Update  Year 

) 

C 

Brother 

1-10 

A10 

Higher  Altitude 

Satel 1 i te  Name 

? 

Sister 

11-21 

A10 

Lower  Altitude 

Satel 1 i te  Name 

Table  5.  Setup  of  Input  Cards  for  Program  FTPRNT 


1 27 


COEFFICIENT  SET  JSEO  IGaMl«»7SI  UPQA  T  EO  TO 


■ 


128 


Figure  C.  Output  Listing  for  Program  FTPRNT 


/ 


4.3.6  Program  Restrictions  and  Timing 

When  the  higher  altitude  satellite  passes  over  the  geomagnetic  north  or 
south  pole,  the  magnetic  field  package  MGFLD2  is  unable  to  trace  the  field 
lines.  Therefore,  an  internal  check  is  made  to  determine  when  this  condition 
exists.  When  it  does,  an  approximation  is  made  for  conjugate  point  locations. 
Whenever  a  conjuncton  is  predicted  using  these  approximations,  a  special 
message  is  written  that  these  results  may  be  unreliable. 

The  lower  the  altitude  of  the  low  altitude  vehicle,  the  more  rapidly  its 
latitude/longitude  coordinates  vary  with  time.  Thus,  the  lower  the  altitude, 
the  smaller  the  time  increment  in  the  ephemeris  calculation  that  is  required  to 
achieve  a  given  accuracy  in  subroutine  SEARCH.  This  translates  into  more  data 
to  be  processed  and,  as  a  result,  increased  central  processor  times  for  lower 
altitude  vehicles. 

4.3.7  Control  Cards  and  System  Requirements 

A  typical  series  of  control  cards  is  as  follows: 


LBOB,  T1000 ,  CM75000. 

ATTACH,  TAPE1 ,  IGRF ,  ID=GEDDES,  MR=1. 
ATTACH,  TAPE2,  SCATHA,  ID=LEWIS,  MR=1. 
ATTACH,  TAPE3,  DMSP  F2,  ID=LEWIS,  MR=1. 
FTN ,  SL,  R=3. 

LDSET ,  PRESET=ZERO. 

LGO ,  PL=20000. 


A  typical  FTPRNT  run  consumes  1000  seconds  of  central  processor  time  and 
requires  65K  octal  core  memory. 
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The  conducting  path  formed  by  the  magnetic  field  line  intersecting  the  SCATHA 
vehicle  is  capable  of  conducting  currents  down  to  each  of  the  conjugate 
footprints  in  the  lower  ionosphere.  Here,  if  ionization  produced  by  solar 
radiation  is  present,  the  conducting  channel  can  extend  horizontally  in  the 
lower  ionosphere.  It  is  of  interest  to  examine  the  solar  illumination  con¬ 
dition  at  these  footprints  to  determine  whether  this  ionospheric  conduction 
mechanism  is  available.  FOOT  has  been  developed  to  make  this  determination. 

4.4.2  Approach 

Figure  6  illustrates  the  flow  of  operations  for  program  FOOT.  The  satellite 
ephemeris  is  accepted  in  LOKANGL  TAPE3  format.  Determination  of  footprint 
location  can  be  made  using  either  of  two  magnetic  field  packages:  MGFLD2  or 

BFLD.  The  latter  uses  the  Olson-Pfitzer  magnetic  field  model.  The  advantage 

claimed  for  this  model  is  that  the  magnetic  field  is  represented  more  realis¬ 
tically  at  the  higher  altitudes  where  currents  due  to  particle  flow  (rather 
than  Earth  core  effects)  predominantly  determine  the  field. 

With  the  latitude,  longitude,  and  altitude  of  the  conjugate  footprints  avail¬ 
able,  Subroutine  SOLLUN  is  called  to  evaluate  the  solar  elevation  angle  at 

these  locations.  A  simple  correction  is  then  made  to  account  for  refraction. 

4.4.3  Input/Output 


The  input  card  setup  for  FOOT  is  shown  in  Table  6.  Table  7  illustrates  a 
typical  set  of  control  cards  for  this  program.  A  sample  output  listing  is 
shown  in  Table  8. 
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pannil  fully  legible  lepioductioa 


Vari  abl  e 
10PT 


COLF  (21 


F  ornat  Description 

13  Selects  Mag  Field 

Model 

0  =  MGFLD2 
1  =  ?L  SPN-PF  I TZF 

2A10  COEF  Label 


F10. 2  Update  Year 


able  (•.  F ornat  of  Input  Cards  for  Program  FOOT 


OF  ERG ,  T 100 ,  CK120000. 

.ATTACH,  TAPE 2,  FOOT,  ID  =  OBERG 
ATTACH,  XX,  IORLIB,  ID  =  WEBER,  FR  =  1. 
LIBRARY,  XX. 

ATTACH,  TAPE  1 ,  IGRF75,  ID  =  WEBER,  MR  =  1 . 
FTN . 

LDSET ,  PRESET  =  ZERO. 

LGO. 


Table  ).  Typical  Sequence  of  Control 
Cards  for  Program  FOOT 


1  't  ;  Z !  >  "i  i  I  •  .•  -  •  *  i  -  -  •  C  -  Z  i  J  t  /  L  t  / 3 1 


.  ^  0  -il  13  3  f  I*)  M  \J  'M  -J  (M  i\J  <'l  <\l  (\l  t\J  f\l  r*-\  T  T  H  U>  X 

>■  J  -I  XJ  ‘  g  4J  J-X  O  <  4  o  ■>  J  o  j  y  o  ;  ’  >  •  4  J  •  J  T  .>  j  -o  >— >  -r 

»  c»  *■*  »»  4j  o  h  "i  x  o  3  **j  X  **•  J  ■  j  .n  n  ;»  rw  r  r*.  n  t-<  r  o  x  •-« 


i\j  »n  x  O  r-  n  j 


„  ^  >u  .\j  N  !*J  w  m  »i  j  ^  j  r  in  ja  in  liN  >iJ  >d  o  i  o_T  J'JI  J>  >j  a  a  a  HTiHHiNjfvjrjr^ 

(N,  1\J  N  ^  ly  (\J  N  f  l  fj  (\J  I'J  rj  (\J  -i  I  j  (\I  ^  "NJ  r\j  (\j  tvj  |\I  rj  N  fj  rg  (\J  ('j  rj  (kj  ^  I  J  W  n  W)  r,  V)  ^  n  »)  ^  M 

Ji  »  'J  n  W  >«  '■'  i;  I J  ui  o  vj  H.  I J  iJ  -J  J-  o  — >  >  o‘  >-l  J*  •>  f»  N *  J  J  r>  i  .>1  r*;  u»  J  ^  'I  vJ  «  S  -O  •‘H  >4  u  41  o*  <  4  r-f 
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3  J  «-  un  rv»  a.  .  j?  i  x  -3  r  3'  U'  ~-J  x  • 

\  .1  1J  T  k|  J  kj  >4  J  .»N  J  'J  'l>  >V 


..  J  j  u»  t~i  '■*•  on  o  «>  *-*  -o  «  j  _>»  o  -4  kkl  -»  *4  »J  I  '"Jy-l  _J  J<  O  w  ,ii  1  J  f>  '\i  ri  j  1_<  •,  >  o  O  ,v-  -O  O  *•»-•*  J  r  -o  *i  <•  «'j 

•  i-j  *v  fj  r*  *-«  •  3  •  3  U‘  -J»  .J»  O  *-  CC  p«*  r^.  Jo  0  -0  o  -O*  1*  «.*<-*  U*  •'  a'  O  •  a>  J*i  U '  T  J  T  J  I  X  J  -T  T  J  f  J  J  X  S 

I  r>;  f  -  O  !•>  /•»  fO  «'J  >\l  .NJ  i\J  >J  ■'  j  k\4  \J  1NJ  I'sl  ‘M  (\J  IVj  k\J  M  "J  l<J  "O  IVi  k\>  -  J  M  '  -  r  J  W  i  \J  fj  «'J  VJ  C-J  I  >  i  >a  '\<  ^  l  \»  4  '\i  '  ’  J  \<  *  si  w 

fs.  xi  j  f  n  ro  "i  n  O'  ii  <*i  r  j  h  h  ('j  'j*  a  *j  n  r  r  o  u>  3  i\<  >v/  a  -j  o  t  y  fv»  r  j  ^  ^  x  t  r  h  o*  i  'vj 

Ok'4f>l*)'*r'.tVJu«\0  t’  U»  ’•*  1J"  U\  J  ~S  •  J*  I  J  ■«!  .1  J  kll  dl  ul  U'  J-f.kJ  J  .  '  ■>'  ->  >  ”  I  *“k  -1  O  .?  i/  A.'  •* »  3  »■*  O  J  »■'  -O 

-O  O  o  »».  vu*  f\  f  *• )  1 SJ  — I  j'  O'*-1  T  M  H  "  U'  '-)  M  uA  /»j  .-i  rr\  u  j  «\j  -%  a  "J  W  O'  J  T  4  0^  yt  *>  O  O  «  'O 


O  L'  O  D  C  O  O  L>  O  r  t!  O  c  '<>  »'  O  ^  IJN  t  T  T  r  J1  *0  -O  rO  ‘ 


-U  -Ni  •»  <r>  Ik.  J  r*  « 


/  r  r»>  o  f 

;'3*o- 

O  Ok.  W  4  l 


t!  O  O  X>  X)  U  o  -k 


J  ^  '1  1*1  lJ>  Uk  t\J  i!>  f  ">  *)  «.'  U\  X  L  >  O  U>  *1  r>  'o  ''I  U'  W'  J  a  H  O  >t)  I|>  (*>  «J  « 

•  •«  r»*  ^  t-  ,1  o  ^  ‘j  i.  ^  4  tk,  k  r  t*  i  t'  wi  f>  t  k.  n  .\  ko  ko  n  «-j  .k  J  r-  i»  > 

ik  j  tj  _»  j*  o  u\  J  >\j  H  o  -y  p»*  j;  y\  f*>  i*i  »-«  oa  j’  s.  yk  j  «•;  w*  .a  o'  o  a>  -J-  * 


J  JJJ’O  ro 

a  ocjoacxTjr— 

«kj  -4  'kj  >g  «\J  *4  .4  '4  4  t.  rj 


i*/  J  u»  u-  --3  k*4  ij  >1  f»  t.M  l»>  ip  T  1 

/  r  T  O  O  »-«  ■>  ■->  k.  1  o-  f 

•"O'3  *v.  uk  I  f;  '4  o  O  ^  jJ  I 


4.1  J  o  J  O  ^  J 


L.-.  -\  J r  3  3  t  t 


*)  4.  O  l»1  3  J  J  i 


r-i  .  j  _  o 

ft:  ^  4-  r^-  'k-  K. 

I  ■  4  i  j  J  l\J  «  k.  f  * 


i  r*»  r»»  '.\j  tk  .j  ««•  ip  u  u  r»^  3  iy  fj  i  >.  t\i  j  o  o  •*<  'Ai  rj  iv  ri  r  u\  .V  o  i>  ^  o  u 

i  is  N  j  M  C  > '  J 1  J  J  /  j  »h  /r^cx  _  o  -o  j  ^  »*•)  *.  J>  cj  _j  -k»  r  _•  w'  i 

U  0.1  I  »«■  .  -  ♦  r*  k'-  o  (  J  I  I*'  o  ,«  Jk  «-)  •«  •■<  4  '  •'-<  *  ‘  *>  ■'  3  J  IN.  ■  J  «  1  I 

.  «v.  o  »  i»  L.  J  N  •'  C  »i  -3  0  4  J  .4  —  -  4  -  -'  *«  O  -  .  >  f*.  m  J*  k.  «-•  J 

.  .  k  :•.  >  u'  '  \:  T  J-  X  t  *  t  *-  -.1  ^  fk  ,  f.  — 

,  ■  i  ;v  »  •  ~J  fj  i’.»  \  IJM  f.  rj  *J  J  4  i  c  J  -  I  -  1  >  I  '  J  J  r  4  j  ■  .  r.  . , 


S  o  jk  J  J  v  >  ,J\  J 


.*•  r ,  rv  'v  •-  .' 


i  fj  n  rs.  •-  t  .r  . 


k-'  r-~  r--  fk  j  r-  ro  r  k-  ^  rj  f  i  OJ  %  ’j  r  i  i\j  r-  .4  r>  '■)  C3  Tk  ^  ~~  C  ■ .  t  X  «•  n.  N.  Is-  *s.  o 

t  t.  n  *■»  •>  i'  o  o  j  *'  "i  er  «.  -it  x  :  «>  *1  1  »■<  ■  0  o  ">  T  '  ^  r*-  ^  ^  ^  N-  ^  ^  >>» 
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•o  00  r  o  o  r»  c  c  -r  -i  oca'  c  o  a  -r  o  >r  a  c  >r  <.1'  o  a*  o  ^  c  kO  x  -O  x  a-  c  x  ^o  u.)  0  kW'  r  .c  x 
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Copt  available  to  DTIC  does  not 
permit  fullY  legible  reproduction 


Table  8.  Printout  of  Solar  Elevation  Angle  at  Footprints  of  S 3- 2  Satellite 


4.5 


Astro-Geophysical  Data  Display 
4.5.1  Introduction 


A  variety  of  astro-geophysical  (AG)  data  consisting  of  the  time  variation  of 
parameters  which  depict  the  state  of  the  magnetosphere  are  collected  by  the  Air 
Force  Global  Weather  Center  (GWC)  and  forwarded  to  AFGL  for  geophysical  analy¬ 
ses.  The  data  are  assembled  onto  tape  at  GWC  using  a  36-bit  word  UNIVAC 
1110  computer.  The  data  are  unpacked,  converted  to  60-bit  words  for  the 
CDC  6600,  and  displayed  on  microfiche. 

4.5.2  AG  Data  Tape  Convention 

The  AG  data  is  written  on  a  7-track  tape  by  a  UNIVAC  computer.  The  following 
definitions  apply:  36  bits  =  12  octal  characters  =  6  field  data  characters  =  1 
word.  The  data  is  written  in  card  images  in  field  data  format;  see  Table  9  for 
the  octal  to  character  conversion.  There  are  70  files  on  the  tape  (7  days  X  10 
data  types).  Each  file  ends  with  a  @E0F  word.  Two  control  words  appear  at  the 
beginning  of  the  tape  and  one  at  the  very  end,  following  the  last  @E0F .  In 
addition,  each  card  image  is  preceded  by  a  control  word.  These  card  control 
words  are  of  the  format 


where 


00AAB6100000 

AA  =  octal  number  of  words  in  the  following  card  image 
BB  =  octal  number  of  words  in  the  preceding  card  image 


The  first  card  image  on  the  first  file  of  the  tape  is 

THIS  TAPE  CREATED  ON  MMDDYY  AT  HHMMSS  BY  AFGWC  FOR  SAMTEC  SCATHA  PROJECT 

This  statement  is  followed  by  the  card  image 

DDMMYY 

which  gives  the  date  on  which  the  data  on  the  next  10  files  was  taken. 

Similarly  the  11th,  21st,  31st. ..files  are  also  preceded  by  a  date  word. 

If  the  tape  was  created  on  Sep  9,  the  date  on  the  first  file  will  be  Sep  8,  the 
date  on  the  11th  file  will  be  Sjp  7,  and  so  on. 
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UNIVAC  UNIVAC 


OCTAL 

CHAR 

OCTAL 

CHAR 

00 

@ 

40 

) 

01 

[ 

41 

- 

02 

] 

42 

+ 

03 

# 

43 

< 

04 

A 

44 

= 

05 

SP 

45 

> 

06 

A 

46 

& 

07 

B 

47 

$ 

10 

C 

50 

★ 

11 

D 

51 

( 

12 

E 

52 

% 

13 

F 

53 

: 

14 

G 

54 

15 

H 

55 

| 

16 

I 

56 

» 

17 

J 

57 

/ 

20 

K 

60 

0 

21 

L 

61 

1 

22 

M 

62 

2 

23 

N 

63 

3 

24 

0 

64 

4 

25 

P 

65 

5 

26 

Q 

66 

6 

27 

R 

67 

7 

30 

S 

70 

8 

31 

T 

71 

9 

32 

U 

72 

1 

33 

V 

73 

9 

34 

W 

74 

/ 

35 

X 

75 

• 

36 

Y 

76 

37 

Z 

77 

f 

Table  9.  Octal  to  Character  Conversion  (UNIVAC) 


I 


The  ten  data  types  (files)  are: 

(1)  Astrogeophysical  Data  Base  (AGDB)  code  type  10  (Raw  Magnetometer  Data) 

(2)  AGDB  code  type  11  (3Hr  Ap) 

(3)  AGDB  code  type  22  (AK  Values) 

(4)  AGDB  code  type  20  (Q  index) 

(5)  GOES  magnetometer  data  (Veh  1) 

(6)  GOES  magnetometer  data  (Veh  2) 

(7)  GOES  particle  data  (Veh  1) 

(8)  GOES  particle  data  (Veh  2) 

(9)  CPA  particle  data  (Veh  6) 

(10)  CPA  particle  data  (Veh  4) 

Figure  7  illustrates  the  organization  of  sequential  files  on  the  AG  tapes. 

Data  for  the  most  recent  day  (i.e.,  the  seventh)  comes  first.  This  data  is 
contained  in  ten  data  files,  one  for  each  data  type.  Next  (i.e.,  the  second 
row  of  the  matrix)  come  the  ten  files  for  the  next  to  last  day.  These  files 
are  numbered  eleven  through  twenty,  inclusive.  This  pattern  repeats  until  all 
ten  files  for  the  earliest  day  have  been  covered.  (They  are  numbered  61 
through  70).  As  an  example  of  the  use  of  the  matrix,  observe  that  the  file 
number  for  GOES  particle  data  (Vehicle  2)  for  the  second  day  is  58.  Note  that, 
although  data  for  individual  days  occurs  in  inverse  chronological  order  on  the 
tape,  the  data  (of  a  given  type)  for  each  day  proceeds  in  normal  chronological 
order. 

Two  types  of  tapes  are  produced  by  GWC:  a  type  1  tape  (old  format)  and  a  type 
2  tape  (new  format).  Type  1  tapes  have  variable  image  lengths  which  are 
allowed  to  span  two  physical  records.  Multiple  AG  files  within  records  can 
also  occur  since  an  @E0F  does  not  correspond  to  a  tape  EOR. 

Type  2  tapes  have  uniform  blocking  of  15  images  per  tape  record.  The  first 
image  in  each  record  is  13  36-bit  words  and  is  blank  filled.  The  remaining 
14  images  are  14  36-bit  words.  After  an  @E0F  the  rest  of  that  image  and 
record  are  blank  filled.  Therefore,  the  next  AG  file  starts  with  the  second 
image  of  the  next  record. 
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DATA  TYPE 


Figure  7.  Sequential  File  Matrix 


4.5.3  Convert  UNIVAC  to  CDC  Code 


The  physical  AG  Tape  record  length  is  224  36-bit  words  in  a  packed  configu¬ 
ration.  Each  6-bit  byte  is  unpacked  from  a  bit  stream.  Character  codes  are 
formed  and  examined  for  @E0F.  Thru  a  selection  of  program  routines,  the 
UNIVAC  equivalent  words  are  formed  and  arrays  of  60-bit  CDC  words  (designated 
MW  in  the  software)  are  unpacked.  Table  10  shows  the  bit  locations  for  36  bit 
packed  words  used  in  Tables  11-14.  These  tables  summarize  the  word  definitions, 
MW  array  numbers,  the  6-bit  words,  and  provide  a  reference  to  the  subroutine 
in  Program  (AGDATA)  which  lists  the  data  as  a  function  of  time. 

The  breakdown  of  the  GOES  magnetometer  data  (File  5  +  o)  and  GOES  particle 
data  (File  7+8)  are  shown  in  Tables  15  and  16.  Tiles  5+7  are  for  Vehicle 
1  data  and  files  6+8  for  Vehicle  2  data.  Each  of  these  files  consists  of 
an  array  of  words  for  six  ten  minute  periods  or  1  hour's  cycle  worth  of 
data.  For  file  5  +  6  a  cycle  is  56  words;  file  7  +  8  is  92  words.  The  bit 
packing  for  each  of  these  files  is  identical. 

The  format  and  the  content  for  the  CPA  particle  data  (Files  9  +  10)  differs 
in  that  all  the  data  is  decoded  as  decimal  numbers.  Table  17  is  a  summary  of 
the  energy  range  of  the  57  channels  of  CPA  data. 


HI 

= 

bits 

1-1 

S2  = 

bi  ts 

7-12 

H2 

= 

bits 

19-36 

S3  = 

bits 

13-18 

T1 

= 

bits 

1-12 

S4  = 

bi  ts 

19-24 

T2 

= 

bits 

13-24 

S5  = 

bi  ts 

25-30 

T3 

= 

bits 

25-36 

S6  = 

bi  ts 

31-36 

Q1 

= 

bits 

1-9 

W  = 

bits 

1-36 

Q2 

= 

bits 

10-18 

H  : 

half 

word 

Q3 

= 

bits 

19-27 

T  : 

th  i  rd 

1  word 

Q4 

bi  ts 

28-36 

0  : 

quarter  word 

SI 

= 

bits 

1-6 

S  : 

sixth 

word 

W  : 

full 

word 

10.  Bit  Definitions  for  Packed  Words 


Table 


MW 

Arra 


6  bit 
word  # 


1st  Image 


WORD  1 

WORD  2 

WORD  3 
WORD  4 
WORD  5 

WORD  6 

WORD  7 

WORD  8 

WORD  9 


HI  =  Word  Count  =  9 

1 

1+6 

S4  =  Day 

2 

S5  =  Hour 

3 

S6  =  Min 

4 

HI  =  Code  Type  =  10 

5 

1+18 

H2  =  AFGWC  Station  No. 

6 

W  =  0 

W  =  WMO  Number 

7 

1+42 

Q1  =  Year  (00-99) 

8 

1+54 

Q2  =  Month 

9 

Q3  =  Status  of  Report 

10 

Q4  =  Processing  Status 

11 

HI  =  XMIN 

12 

1+66 

H2  =  YMAX 

13 

HI  =  YM1N 

14 

1+6 

H2  =  YMAX 

15 

HI  =  ZMIN 

16 

1  +  18 

H2  =  ZMAX 

17 

HI  =  MAX  GAMMA 

18 

1+30 

Q3  =  ak  value 

19 

Q4  =  k  value 

20 

Table  li.  CODE  TYPE  10  (Magnetometer  Data)  -  Subroutine  FILE  1 


1 31' 


MW  6  bi t 

Array  word  £ 

WORD  1:  Hi  -  Word  Count  =5  1  1+6 

54  =  Day  2 

55  =  Hour  3 

56  =  Min  A 

WORD  2:  HI  =  Code  Type  =  11  6  1+16 

H2  =  AFGWC  Station  No.  6 

WORD  3:  W  =  0 

WORO  4:  HI  =  (2  digit  year)  *100  +  Month  7  1+42 

H2  =  24  hr  ap  8 

WORD  5:  HI  =  3  hrly  ap  9  1+54 

H2  =  Kp  10,11,12 


Table  yi,  CODE  TYPE  11  (Geomagnetic  Indices)  -  Subroutine  FILE  2 
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MW 

Array 


6  bit 
word  # 


1st  Image 


WORD 

1: 

HI 

= 

Word  Count  =  11 

1 

1+6 

S4 

= 

Day 

2 

S5 

= 

Hour 

3 

$6 

= 

Min 

4 

WORD 

2: 

HI 

Code  Type  =  22 

5 

1+18 

H2 

= 

AFGWC  Station  No. 

6 

WORD 

3: 

W  * 

0 

WORD 

4: 

W  * 

WMO  Number 

7 

1+42 

WORD 

5: 

Qi 

s 

Year 

8 

1+54 

Q2 

= 

Month 

9 

Q3 

= 

AK  Value 

10 

Q4 

Status  =  0 

11 

WORD 

6 

Qi 

= 

K.  value  of  First  3  Hr  Period 

12 

1+66 

Q2 

= 

K  value  of  Second  3  Hr  Period 

13 

Q3 

K  value  of  Third  3  Hr  Period 

14 

Q4 

K  value  of  Fourth  3  Hr  Period 

IS 

WORD 

7 

Ql 

K  value  of  Fifth  3  Hr  Period 

16 

1+6 

Q2 

= 

K  value  of  Sixth  3  Hr  Period 

17 

Q3 

= 

K  value  of  Seventh  3  Hr  Period 

18 

Q4 

= 

K  value  of  Eighth  3  Hr  Period 

19 

WORD 

8 

HI 

= 

Indicator  of  Phenomena 

20 

1+18 

H2 

= 

Time  of  Phenomena 

21 

WORD 

9 

HI 

= 

Time  of  Min  Gamma  Value 

22 

1+30 

H2 

= 

Minimum  Gamma  Value 

23 

WORO 

10 

HI 

Time  of  Instantaneous  H 

24 

1+42 

H2 

= 

Instantaneous  H  (Gammas) 

25 

WORD 

11 

HI 

= 

0 

26 

1+54 

H2 

= 

Big  AK  of  AK 

27 

2nd  Image 


Table  13.  CODE  TYPE  22  -  (Geomagnetic  Indices)  Subroutine  FILE  3 
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MW 

Arra 


6  bit 
word  # 


1st  Image 


WORD  1:  HI  =  Word  Count  =  10  1 

54  =  Day  2 

55  =  Hour  3 

56  =  Min  4 

WORD  2:  HI  =  Code  Type  =20  5 

H2  =  0  6 

WORD  3:  W  =  0 

WORD  4:  HI  =  Satellite  ID  7 

H2  =  Rev  Number  8 

WORD  5:  Equatorward  Latitude  (tenths)  9 

WORD  6:  Equatorward  Longitude  (tenths)  10 

WORD  7:  Equatorward  Corrected  Geomagnetic  11 

Latitude 

WORD  8:  Corrected  Geomagnetic  Local  Time  12 

WORD  9:  Q  Index  (Tenths)  13 

WORD  10:  T1  =  H/J  Data  14 

T2  =  Data  Qual ifier  15 

T3  =  Aurora  Characteristic  16 


1+6 


1+18 


1+42 

1+54 

1+66 

1+6  2nd  Image 

1+18 

1+30 

1+42 


Table  14.  CODE  Type  20  -  (Auroral  Index)  Subroutine  FILE  4 
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Words  1-48  are  data  for  six  ten  minute  periods  from  eight  data  lines  of  the 
GOSMM  code. 

words  i-6  Data  line  containing  5  minute  averaged  magnetic  field  total 

magnitude  values  valid  on  the  hour  and  every  10  minutes  during  the 
following  hour. 

Words  1-11  Data  line  containing  5  minute  averaged  values  of  the  x-component 
of  the  magnetic  field  in  an  earth-centered  solar-ecliptic  rect¬ 
angular  coordinate  system  having  its  x-axis  oointed  toward  the 
sun,  its  z-axis  pointed  northward  perpendicular  to  the  ecliptic 
plane  and  its  y-axis  in  the  ecliptic  plane. 

/» era s  Data  line  containing  similar  5  minute  averages  values  of  the 

/-component  of  the  magnetic  field. 

words  ; 0  -  a 4  Data  line  containing  5  minute  averaged  values  of  the  z  comparer* 
of  the  magnetic  ’■.eld. 

/.<■'  cs  >  -i-.'  Data  line  containing  standard  deviations  of  the  averaged  total 
magnetic  field  magnitude  for  each  report  during  the  near. 

-eras  31-36  Data  line  containing  similar  standard  deviation  values  for  the 
x-component  of  the  magnetic  field. 

words  37-42  Data  line  containing  similar  standard  deviation  values  for  the 
y- component  of  the  magnetic  field. 

Words  43-48  Data  line  containing  similar  standard  deviation  values  for  the 
z-component  of  the  magnetic  field. 

Words  1-48  are  packed  with  the  FLD  function  as  follows: 

0  13  17  27  29  35  Starting  bit  of  field 

Unused  Qual  Field  Sign  Power  of  10  Data  Contents 


Word  49  Local  time  at  sub-satellite  point  on  the  hour. 
Word  50  Latitude  sign  (1  is  for  North,  2  for  South) 

Word  51  Latitude  of  sub-satellite  point,  (whole  degrees) 
Word  52  Longitude  sign  (1  is  for  West,  2  for  East) 

Word  53  Longitude  of  sub-satellite  point,  (whole  degrees) 
Word  54  Year-Month-Day  of  data  for  this  hour  (YYMMDD). 
Word  55  Currently  unused. 

Word  56  Currently  unused. 


Table  15.  GOES  Magnetometer  Data  -  Subroutine  FILE  5 
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Words  1-84  are  data  for  six  ten  minute  periods  from  14  channels  of  the 
GOSPP  coded  data. 


Particle  Type 


Energy  Range 


Words 

1-6 

ELECTRONS 

>  2 

MEV 

Words 

7-12 

PROTONS 

0.8-4 

MEV 

Words 

13-1 8 

PROTONS 

4-8 

MEV 

Words 

19-24 

PROTONS 

8-16 

MEV 

Words 

25-30 

PROTONS 

16-215 

MEV 

Words 

31-36 

PROTONS 

36-215 

MEV 

Words 

37-42 

PROTONS 

80-215 

MEV 

Words 

43-48 

PROTONS 

215-500 

MEV 

Words 

49-54 

ALPHAS 

4-10 

MEV 

Words 

55-60 

ALPHAS 

10-16 

MEV 

Words 

61-66 

ALPHAS 

18-56 

MEV 

Words 

67-72 

ALPHAS 

71-150 

MEV 

Words 

73-78 

ALPHAS 

167-245 

MEV 

Words 

79-84 

ALPHAS 

340-392 

MEV 

Words 

1-84  are  packed  with 

the 

FLD  function  as  follows: 

0 

13 

17 

27 

29 

35 

Starting 

Unused  Qual 

Field 

Sign 

Power 

of  10 

Data  Con 

Word  85  Local  time  on  the  hour. 

Word  86  Latitude  sign  (1  is  North,  2  for  South) 

Word  87  Latitude  in  degrees 

Word  88  Longitude  sign  (1  is  West,  2  is  East) 

Word  89  Longitude  in  degrees 

Word  90  Year-Month-Day  of  data 

Word  91  Unused. 

Word  92  Unused. 


Table  16.  GOES  Particle  Data  -  Subroutine  FILE  7 
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Satel 1 ite 
Look  Ang 


Particle 


+60 

DEG 

ELECTRONS 

30  - 

300  KEV 

+60 

DEG 

ELECTRONS 

44  - 

300  KEV 

+60 

DEG 

ELECTRONS 

64  - 

300  KEV 

+60 

DEG 

ELECTRONS 

95  - 

300  KEV 

+60 

DEG 

ELECTRONS 

139  - 

300  KEV 

+60 

DEG 

ELECTRONS 

204  - 

300  KEV 

+30 

DEG 

ELECTRONS 

30  - 

300  KEV 

SAT 

EQ 

ELECTRONS 

30  - 

300  KEV 

SAT 

EQ 

ELECTRONS 

44  - 

300  KEV 

SAT 

EQ 

ELECTRONS 

64  - 

300  KEV 

SAT 

EQ 

ELECTRONS 

95  - 

300  KEV 

SAT 

EQ 

ELECTRONS 

139  - 

300  KEV 

SAT 

EQ 

ELECTRONS 

204  - 

300  KEV 

-30 

DEG 

ELECTRONS 

30  - 

300  KEV 

-30 

DEG 

ELECTRONS 

44  - 

300  KEV 

-30 

DEG 

ELECTRONS 

64  - 

300  KEV 

-30 

DEG 

ELECTRONS 

95  - 

300  KEV 

-30  DEG 

ELECTRONS 

139  - 

300  KEV 

-30 

DEG 

ELECTRONS 

204  - 

300  KEV 

-60 

DEG 

ELECTRONS 

30  - 

300  KEV 

-60 

DEG 

ELECTRONS 

44  - 

300  KEV 

-60 

DEG 

ELECTRONS 

64  - 

300  KEV 

-60 

DEG 

ELECTRONS 

95  - 

300  KEV 

-60 

DEG 

ELECTRONS 

139  - 

300  KEV 

-60 

DEG 

ELECTRONS 

204  - 

300  KEV 

ELECTRONS 

0.20  - 

2.0  MEV 

ELECTRONS 

0.29  - 

2.0  MEV 

ELECTRONS 

0.43  - 

2.0  MEV 

ELECTRONS 

0.63  - 

2.0  MEV 

ELECTRONS 

0.98  - 

2.0  MEV 

ELECTRONS 

1.36  - 

2.0  MEV 

PROTONS 

50  - 

500  KEV 

PROTONS 

63  - 

500  KEV 

PROTONS 

80  - 

500  KEV 

PROTONS 

100  - 

500  KEV 

PROTONS 

126  - 

500  KEV 

PROTONS 

158  - 

500  KEV 

PROTONS 

198  - 

500  KEV 

PROTONS 

252  - 

500  KEV 

PROTONS 

316  - 

500  KEV 

PROTONS 

396  - 

500  KEV 

PROTONS 

0.3 

0.4 

(0.35) 

MEV 

PROTONS 

0.4 

0.53 

(0.46) 

MEV 

PROTONS 

0.53 

0.71 

(0.62) 

MEV 

PROTONS 

0.71 

0.94 

(0.82) 

MEV 

PROTONS 

0.94 

1.25 

(1.1) 

MEV 

PROTONS 

1.25 

1.66 

(1.46) 

MEV 

PROTONS 

1.66 

2.80 

(1.93) 

MEV 

Table  17. 

CPA  Particle 

Data-Subroutine 

FILE  9 
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Particle 


Type 

Enerqy  Range 

PROTONS 

2.80  - 

4.73 

(3.76)  MEV 

PROTONS 

4.73  - 

8.00 

(6.36)  MEV 

PROTONS 

8.00  - 

13.5 

(10.8)  MEV 

PROTONS 

13.5  - 

22.8 

(18.2)  MEV 

PROTONS 

22.8  - 

33.2 

(28.0)  MEV 

PROTONS 

33.2  - 

48.4 

(40.8)  MEV 

PROTONS 

48.4  - 

70.6 

(59.5)  MEV 

PROTONS 

70.6  - 

103. 

(86.8)  MEV 

PROTONS 

103. 

150. 

(126.)  MEV 

Table  17. 

CPA  Particle  Data- 

Subroutine 

FILE  9  (continued) 
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Five  Fortran  programs  process  AG  related  information.  Table  18  is  a  list  of 
these  programs  along  with  a  summary  description  of  each. 


Run  Time 


PROGRAM 

Core 

(Sec) 

DESCRIPTION 

LOCZD 

60K 

500 

Dump  AG  Tape  (complete) 

LOCZL 

60K 

2 

List  AGMMDD  Tapes 

LOCZR 

60K 

500 

Reformat  AG  Type  1  (old) 
create  AG  Type  2  (new) 

tape. 

LOCZS 

135K 

1800 

List  data  and  produce  microfiche 
plots  from  AG  type  2  (new)  tape 

LOCZT 

60K 

2 

Dump  AG  Tape  (1st  record)  header 

TABLE  18:  AG  Programs 

Program  LOCZS  with  its  system  of  subroutines  converts  AGMMDD  type  2  data 
tapes  from  UNIVAC  formatted  code  to  CDC  code.  From  an  input  punched  card,  the 
number  of  days,  the  total  number  of  files  utilized,  and  file  selection  order 
is  determined.  An  example  of  a  file  1  chronological  selection  would  be:  71, 

61,  51,  41,  31,  21,  11,  1.  As  the  AGMMDD  tape  is  processed  @E0F  detection 
defines  the  file  partitions.  Thus  all  AG  data  and  file  markers  with  the  excep¬ 
tion  of  the  CPA  data  is  contained  in  a  single  file.  The  CPA  data  is  placed  in 
14  addressable  files.  By  volume  84%  of  the  AG  data  is  CPA  particle  data.  This 
file  structure  system  permits  a  considerable  run  time  savings  since  the  search 
for  data  is  greatly  reduced.  Figure  8  illustrates  the  flow  of  information  for 
Program  LOCZS. 

Through  a  system  of  frames  which  contain  up  to  four  plots  each,  190  plots  are 
contained  in  53  frames.  A  summary  of  these  plots  is  shown  in  Table  19.  AG 
data  files  of  the  same  type  must  be  chosen  sequentially  with  no  restriction 
as  to  chronological  ordering.  If  a  file  is  empty  the  day  counter  will  not  be 
activated  because  the  AG  tape  contains  70  files,  full  or  not.  Since  the  day 
counter  actually  checks  EOF  status,  a  frame  will  be  produced  only  if  data  is 
available  to  be  plotted. 


FROM  PAGE  ! 


3  ■  4  PIOTS/FRAME 
1  -  2  PIOTS/FRAHE 

4  FRAMES  TOTAL 


Figure  8b.  Flow  Chart  for  LOCZS  (Continued) 
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AG  File  # 

F  rames 

Plots/Frame 

Description 

1 

7 

4 

Raw  Magnetometer  Data 

2 

1 

1 

3  Hr  Ap 

4 

1 

1 

Q  Index 

7 

3 

4 

GOES  Particle  Data  (Vehicle  1 

1 

2 

8 

3 

4 

GOES  Particle  Data  (Vehicle  2 

1 

2 

9 

14 

4 

CPA  Particle  Data  (Vehicle  6) 

1 

1 

2 

4 

Angular  Distribution  and 

High /Low  Energy  Summaries 

1 

1 

10 

14 

4 

CPA  Particle  Data  (Vehicle  4) 

1 

1 

2 

4 

Angular  Distribution  and 
High/Low  Energy  Summaries 

1 

1 

TABLE  19.  AG  Plot  System 

For  most  types  of  data,  the  values  plotted  are  essentially  the  raw  numbers 
unpacked  from  the  data  tape.  Magnetometer  data,  though,  is  an  exception. 

Table  11  shows  that  both  maximum  and  minimum  values  are  reported  for  each  of 
the  X,  Y,  and  Z  coordinates.  This  data  occurs  at  the  rate  of  one  set  of  values 
every  90  minutes.  The  value  plotted  at  any  given  time  is  derived  from  both  the 
value  quoted  at  that  time  and  the  value  corresponding  to  the  time  90  minutes 
earlier.  The  minimum  is  taken  for  the  minima  at  these  two  times.  Similarly, 
the  maximum  is  taken  for  the  two  maxima.  The  number  plotted  is  obtained  by 
subtracting  the  resulting  minimum  of  minima  from  the  corresponding  maximum  of 
maxima.  This  procedure  is  followed  for  each  of  the  three  rectangular  compo¬ 
nents  of  the  field. 
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Programs  LOCZD  and  LOCZT  are  similar  in  that  each  examines  an  AG  tape  as  to 
type  and  creation  date.  LOCZD  has  additional  logic  that  checks  an  entire 
tape  for  overall  file  content,  and  has  an  option  for  an  octal  dump.  LOCZT, 
on  the  other  hand,  checks  only  the  first  record  for  type  ,  creation  date,  and 
file  information  to  delineate  the  AGMMDD  tape  label,  where  MM  is  the  month 
and  DD  is  the  day  of  the  first  chronological  file. 

Program  LOCZL  lists  the  AGMMDD  tapes  at  AFGL  along  with  their  inclusive  file 
dates  and  types. 

Program  LOCZR  creates  an  AG  type  2  (new)  tape  from  an  AG  type  1  (old)  tape. 
Typically,  type  2  tapes  are  received  at  AFGL  weekly,  despite  the  fact  that 
type  1  processing  is  available.  A  parallel  programming  task  was  abandoned 
when  it  was  decided  not  to  continue  type  1  reduction.  It  was  determined  that 
the  type  2  tape  format  was  more  conducive  to  efficient  chronological  sorting. 
However,  when  only  the  type  1  tape  was  received,  and  since  GWC  retains  AG  data 
for  only  60  days,  program  LOCZR  provides  a  means  to  reduce  significantly  more 
data  than  would  otherwise  be  possible. 


152 


4.6  References 


1.  Hess,  W. 
1968. 

2.  Alfven, 
Oxford  Univ. 


N.  "The  Radiation  Belt  and  Magnetosphere,"  Blaisdell,  Waltham,  MA, 

H.,  and  C.-6.  Falthammar,  "Cosmical  Electrodynamics , "  2nd  Ed, 
Press,  London,  1963. 


153 


5.0  Rocket  Trajectory  System 

This  section  discusses  updates  made  to  the  rocket  trajectory  system  since 
preparation  of  a  comprehensive  report  and  user's  guide.*  These  changes 
include: 

(1)  Addition  of  dynamics  models  for  new  vehicles; 

(2)  Upgrading  of  Ft. Churchill  boresight  directional  corrections  and 
application  of  same  to  other  radars; 

(3)  Conversion  of  White  Sands  trajectory  reports  to  AFGL  TAPE4  format, 

b.l  Rocket  Dynamics  Models 

Thrust,  mass  and  drag  characteri sties  of  several  new  vehicles  have  been  added 
to  DRIVEB  since  publication  of  Ref.  1.  These  are  listed  in  Table  1.  Test 
integration  and  filtering  runs  resulted  in  adjustments  in  some  cases. 

5.1.1  Multiple  Modules  and  Thrust  Angles 

In  some  cases  the  vehicle  separates  into  modules,  which  are  then  to  be  analyzed 
individually  by  DRIVEB.  Furthermore  the  separate  module(s)  may  undergo  attitude 
maneuvers  resulting  in  thrusting  at  angles  other  than  the  velocity  direction. 

An  example  of  same  is  the  ARIES  vehicle  from  which  the  payload  separates  63  sec 
after  launch.  At  90  sec,  there  is  further  separation  of  payload  into  sensor 
and  target,  the  latter  of  which  thrusts  at  various  angles  off  the  velocity 
vector  during  the  remainder  of  the  flight.  Furthermore  the  boost  phase  is 
known  to  iiave  non-zero  velocity-thrust  angles.  In  this  case  the  thrust  direction 
has  been  defined  by  the  angle  from  vertical.  It  is  assumed  to  lie  in  the 
azimuthal  plane  of  the  flight. 

The  procedures  for  handling  these  are  are  follows:  prior  to  separation,  sensor 
and  target  modules  are  assigned  the  booster's  dynamic  model.  Thus  separate 
information  for  these  models  has  been  added  only  after  separation. 
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Rocket  _Key 

Sergeant-Hydac  02 

Nike-Orion  12 

Taurus-Orion  13 

Aries  Sensor  26 

Aries  Target  27 

Talos  -  Castor  51 


Table  1.  Dynamics  Models  for  New  Rockets 
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A  special  data  array  has  been  added  to  hold  thrust  angle  as  a  piecewise 
linear  function  of  time.  This  is  applied  to  the  booster  initially  and  later 
to  the  target. 

5. 1.1.1  Suggestions  for  Improved  Handling  of  Multiple  Module  Missions 

The  above  procedure  could  result  in  a  proliferation  of  models  as  m-.1  iple- 
module  missions  become  more  frequent.  Furthermore  separations  cou'  ■.  occur 
with  finite  velocity.  One  way  to  handle  this  latter  occurrence  is  to  build  a 
short  impulse  into  the  thrust  model  (say,  1/2  sec  in  duration).  This  pro¬ 
cedure  has  been  employed  in  some  cases  not  permanently  built  into  tne  program. 
However  this  could  use  up  valuable  storage  space  in  the  model  data  array.  More 
systematic  procedures  are  suggested  as  follows. 

Selective  Tracking  of  Modules  (Payloads  or  Booster)  in  DRIVEB. 

1.  Up  to  5  modules  including  booster. 

2.  Specifications  for  each  module. 

a)  Number  (1-5) 

b)  Weight  (with  and  without  fuel),  area,  drag 

c)  Alphanumeric  identifier  for  thrust  model  (if  any) 

3.  Separation  specifications 

a)  Time 

b)  Contents  of  Vehicle  #1  (module  numbers) 

c)  Contents  of  vehicle  #2  (module  numbers) 

d)  Separation  velocity  (vehicle  #1  -  vehicle  #2) 

Vehicle  #1  will  be  the  front  vehicle. 

Vehicle  #2  will  be  the  back  vehicle. 

4.  Thrust  and  angles  for  payload  modules. 
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Logic: 


1.  For  the  selected  module  the  following  is  saved  after  completion  of 
each  integraton  step. 

a)  Separation  times  immediatedly  before  and  after  current  time,  for 
separations  involving  the  module  (TSEP1  and  TSEP2) 

b)  Modules  to  which  it  is  connected  Mod  (I),  1=1,  NMOD 

c)  Total  Weight,  drag  and  area  of  package  and  thrust  model  identifier 

2.  At  each  new  step  test  to  see  if  time  is  still  between  TSEP1  and  TSEP2; 
if  so  compute  forces  and  perform  integration  step  in  usual  fashion;  if 
not,  update  module  information  mentioned  in  step  1  and  make  velocity 
correction  in  accordance  with  specified  separation  velocity;  then  proceed 
with  usual  force  computation. 

5.2  Bores ight  Corrections 

Corrections  have  been  incorporated  to  account  for  directional  differences 
between  the  RF  beam  and  the  mechanical  axis  of  a  radar.2  The  raw  (azimuth  and 
elevation)  angles  represent  the  mechanical  axis  direction,  while  the  tracked 
target  lies  along  the  RF  beam.  Previously,  corrections  had  been  made  to  Ft. 
Churchill  data  only.  These  corrections  were  independent  of  the  elevation  angle 
of  the  radar  and  hence  were  accurate  for  low  elevations  only. 

To  determine  the  azimuth  and  elevation  corrections  in  the  general  case  we  first 
define  the  unit  vector  components  of  the  RF  direction  in  the  radar-fixed  co¬ 
ordinate  system: 

x  =  sin  K  cos  P 
y  =  cos  K  cos  P 
z  =  sin  P 
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where  K  and  P  are  the  azimuth  and  elevation  of  the  RF  beam  when  the  mechanical 

axis  reads  azimuth=elevation=0.  Thus  the  radar-fixed  y  axis  lies  along  the 

mechanical  axis,  the  z  axis  points  vertically  upward  when  the  mechanical 

elevation  =  0,  and  x-y-z  defines  a  right-handed  coordinate  system.  The  angles  K 

3 

and  P  can  be  found  from  pi unge-and-rotate  methods.  To  obtain  the  RF 
azimuth  and  elevation  in  the  general  case  we  require  the  coordinate  transfor¬ 
mation  from  the  radar-fixed  system  to  the  local-vertical  system: 


where  A  and  E  are  the  mechanical  axis  azimuth  and  elevation  angles.  The  local 
vertical  coordinates  x‘,  y',  z'  of  the  RF  beam  are  related  to  the  RF  azimuth  and 
elevation  angles  A'  and  E'  by: 

x'  =  sinA1  cosE' 
y1  =  cosA'  cosE1 
z'  =  sin  E' 

Carrying  out  the  required  matrix  multiplication  and  substituting  for  x,  y,  z 
from  above  results  in 

cosE'  sinA1  =  cosA  sinK  cos  P  +  sinA  (cosK  cosE  cos  P- 
sinE  sin  P) 

cosE*  cosA'  =  -sinK  cos  P  sinA 

+cosA  (cosE  cosK  cos  P  -  sin  E  sin  P) 

sinE'  =  sinE  cosK  cosP  +  cosE  sinP 

Using  the  trigonometric  identities  for  the  sine  and  cosine  of  the  sum  of  two 
angles  it  is  found  that 


r 


sinK  cos  P 

tan  t  =  _ 

cosKcos  (E  +  P)  -  (1-cosK)  sinE  sin  P 

sinE*  =  cosK  sin  (E  +  P)  +  cosE  sin  P  (1-cosK) 

Neglecting  terms  of  3rd  order  and  higher  in  the  angles  K  and  P  simplifies 
these  expressions  to 


tanK 

tan  t  =  - 

cos(E+  P) 

sin  E*  =  cos  K  sin  (E  +  P) 

These  results  are  exact  if  either  K  or  P  is  equal  to  zero,  the  latter  of 
which  is  often  the  case.  They  have  therefore  been  incorporated  into  DRIVEA. 

It  is  therefore  necessary  to  read  DRIVEA  input  cards  5-7^  regardless  of 
radar  range  parameter  ICODE,  which  is  now  used  only  to  distinguish  Eglin 
(ICODE  =  1)  from  the  others  (ICODE  =  2),  for  the  purpose  of  adding  Eglin's 
flight  line  direction  to  the  Az  angle. 

5.3  WSMR  Trajectory  Report  Conversion 

White  Sands  Missile  Range  (WSMR)  often  makes  available  results  of  its  radar 
data  reduction  on  UNI  VAC  binary  and  BCD  tapes.4  It  is  often  useful  to 
compare  these  results  with  those  obtained  from  the  AFGL  rocket  trajectory 
system.  To  facilitate  this  procedure,  software  was  developed  for  conversion 
of  the  WSMR  3x80  binary  data  to  the  AFGL  rocket  trajectory  report  (TAPE  4) 
format.5  This  software  consists  of  two  programs: 

1.  Program  WSMRBIN,  which  creates  an  intermediate  file  containing 
geocentric  position,  velocity,  and  acceleration  from  a  WSMR  3x80  binary 
tape; 

2.  Program  FORGEN,  which  converts  the  file  written  by  WSMRBIN  to  an  AFGL 
rocket  trajectory  report  tape,  using  an  earlier  format.6  Raw  data,  resi¬ 
duals  and  covariances  are  omitted  (zero-filled). 
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Program  NEWFOR,  which  replaces  FORGEN,  provides  these  variables  in  the  new 
format  given  in  Ref.  1.  The  variables  available  are  indicated  in  Table  2 
(remaining  words  zero-filled).  The  input  tape  (TAPE2)  as  created  by  WSMRBIN) 
consists  of  the  following: 

Header: 

HEAD  (3-4)  launcher  latitude,  longitude  (W)  (deg) 

HEAD  (5-6)  radar  latitude,  longitude(W)  (deg) 

IHEAD  (8-10)  launch  month,  day,  year 
IHEAD  (11)  seconds  UT  launch  time 

Data  Records: 

GMT(sec)  followed  by  geocentric  position,  velocity,  and  accleration 
vectors  (km  , Km/sec,  Km/sec^),  followed  by  altitude  (Km). 

The  following  card  inputs  are  required: 

Card  1  Cols  1-70(7A10)  TITLE (1-7 ) 

Card  2  Cols  1-20 (2A10)  DAT(5-6):  Rocket  id 

TITLE (8)  is  left  blank.  The  radar  and  launcher  altitudes  are  set  to  1.2312 
Km. 


Program  WSMRBIN  is  a  modification  of  program  WSMR?,  an  SUN  library  program 
which  was  written  for  conversion  of  WSMR's  multi-radar  3x80  UN1VAC  1108  binary 
tapes.  The  input  file  (TAPE  13)  is  as  described  in  Ref.  4,  with  the  following 
provisos: 

One  header  record  per  station  plus  one  for  multi-station  solution  if  it 
exists.  One  80  word  data  group  per  station  per  time  line  (+1  for  multi-station 
solution  if  it  exists),  blocked  in  3  groups  per  3x80  physical  record. 

Four  input  cards  are  read,  as  described  in  Table  3.  One  binary  output  file  is 
generated  per  station  (plus  one  for  multi-station  solution)  on  tapes  2  thru  7 
in  the  format  previously  discussed  for  input  to  NEWFOR. 
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Header  record 
Word  # 


Symbol 


Description 


1-8 

TITLE (1-8 ) 

9-10 

DAT (1-2 ) 

18 

DAT (4) 

19-20 

DAT(5-6) 

21-23 

PROG ID ( 1 - 3 ) 

24 

HRAD 

25-27 

RLODEG,  RLOMIN 
RLOSEC 

28-30 

RLADEG,  RLAMN, 

RLASEC 

31 

AALT 

32-34 

ALODEG,  ALOMIN 
ALOSEC 

35-37 

ALADEG,  ALAMIN 
ALASEC 

38 

TLNCH 

43 

DT 

44 

TLI 

45 

TLF 

46 

TSKIP 

80  character  alphanumeric  id 
Run  date,  blank 
Launch  date 
Rocket  id 

Program  id  for  plotted  output 
Radar  altitude  (Km) 

Radar  longitude  (deg  W,  min,  sec) 

Radar  latitude  (deg  N,  min,  sec) 

Launcher  altitude  (Km) 

Launcher  longitude  (deg  W,  min,  sec) 

Launcher  latitude  (deg  N,  min,  sec) 

Launch  time  UT  (sec) 

Computation  time  step  size  (sec) 
Initial  print  time  from  launch 
Final  print  time  from  launch) 

Tape  written  every  TSKIP  secs 


Table  2.  WSMR  Trajectory  Report  (AFGL  Format) 


J 


\ 


J 


I 


L. 
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Data  records 


Word  # 

Symbol 

Description 

1 

GMT 

Universal  time  (sec) 

2 

TIME 

Tine  after  launch  (sec) 

4 

RAG 

Right  ascension  of  Greenwich  (rad) 

5-7 

P  2  ( 1  -  3 ) 

Filtered  geocentric  pos.  vector  (km) 

3-10 

P2 (4-6 ) 

Filtered  geoc.  vel •  vector  (Km/sec) 

11-13 

P2( 7-9) 

Filtered  geoc.  acc.  vector  (Km/sec^) 

14-16 

PVL (1-3) 

Filtered  launcher  ref.  pos.  (Km) 

17-19 

PVL (4-6 ) 

Filtered  launcher  ref.  vel.  (Km/sec^) 

20-22 

PVL (7-9 ) 

Filtered  launcher  ref.  acc.  (Km/sec^) 

23-25 

OVL ( 1-3 ) 

Filtered  launcher  ref.  range,  az,  el  (Km, deg) 

26-28 

OVL (4-6 ) 

Filtered  launcher  ref.  range,  az,  el  rates 
(Km/sec,  deg/sec) 

65-67 

OVR (1-3) 

Same  as  wds.  23-25,  but  for  radar 

68-70 

OVR (4-6 ) 

Same  as  wds  26-28,  but  for  radar 

71-73 

PVR ( 1-3 ) 

Same  as  wds  14-16,  but  for  radar 

74-76 

GPV(l-3) 

Geodetic  alt.,  long(W),  lat.  (Km,  deg) 

77 

VR 

Local  velocity  magnitude  (Km/sec) 

78,79 

AZR ,  ELR 

Azimuth,  elevation,  of  local  velocity  vector  (deg) 

90 

GRANGE 

Ground  range  along  spheriod  (Km) 

91 

ACL MAG 

Launcher  ref.  accel .  magnitude  (Km/sec^) 

Table  2.  WSMR  Trajectory  Report  (AFGL  Format)  Continued 
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Card 

No. 

Variable 

Name 

Card 

Col . 

Format 

Variable  Description 

1 

IRKNO 

1-10 

A10 

ROCKET  NUMBER 

1 

ISTNO 

11-15 

15 

NUMBER  OF  RADAR  STATIONS 

1 

STPTIM 

16-29 

F10.3 

TIME  (SEC  GMT)  to  STOP  PROCESSING 

2 

LMON 

1-2 

12 

MONTH  OF  LAUNCH 

2 

LDAY 

4-5 

12 

DAY  OF  LAUNCH 

2 

LYR 

7-8 

12 

YEAR  OF  LAUNCH 

2 

I  HR 

10-19 

F10.4 

HOUR  OF  LAUNCH 

2 

I  MI  N 

20-29 

F10.4 

MINUTE  OF  LAUNCH 

2 

I  SEC 

30-39 

F10.4 

SECOND  OF  LAUNCH 

2 

IOPT 

40 

11 

1  «  NO  LISTING  -  TAPE  OUTPUT 

2  =  LISTING  -  NO  TAPE  OUTPUT 

3  =  LISTING  -  TAPE  OUTPUT 

2 

IMULT 

41-42 

12 

1  =  MULTI-STATION  SOLUTION  ON  TAPE 

ANY  OTHER  VALUE  =  NO  MULTI-STATION 
ON  TAPE 

3 

IROCNO 

1-30 

615 

RADAR  STATION  ID's  UP  TO  SIX 

INCLUDING  999  FOR  MULTI-STATION 
SOLUTION  IF  PRESENT. 

4 

10(1)- 

10(40) 

1-80 

4012 

WORD  NUMBERS,  IN  80  WORD  DATA  GROUPS, 
OF  VARIABLES  TO  BE  LISTED  IN  PRINT¬ 
OUT  (MAXIMUM  OF  40  ALLOWED) 

Table  3.  WSMRBIN  Punched  Card  Input 
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